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Introduction 


This  project  aimed  at  developing  new  reconstruction  methods  that 
improves  the  quality  and  quantification  of  functional  properties  in 
near  in  near  infrared  (NIR)  imaging  [1-3].  This  NIR  imaging 
technique  uses  NIR  light  (600  nm  -  1000  nm)  to  probe  the  breast 
tissue  delivered  through  fiber  optic  bundles,  which  are  also  used  to 
collect  the  propagated  diffused  light  (fiber-optic  setup  is  shown  in 
Figure  1)  [2,3],  giving  its  name  diffuse  optical  tomography  (DOT). 
Modeling  of  this  diffuse  light  propagation  requires  advanced 
computational  methods  due  to  the  non-linear  nature  of  light 
propagation,  mainly  due  to  the  dominance  of  light  scattering  over 
attenuation  [1],  Most  of  the  current  computational  models  are  in 
two-dimensions  (2D),  even  though  light  travels  in  three-dimensions 
(3D),  due  to  computational  complexity  and  limited  number  of  measurements.  Estimation  of  functional 
properties  is  highly  dependent  on  these  computational  models  and  development  of  new  computational  methods 
(image  reconstruction  methods)  in  3D  that  reduce  computational  complexity  is  the  main  objective  of  this 
project. 

Specifically,  the  aims  of  this  project  are 

1)  Reducing  the  computational  complexity  of  3D  optical  imaging  by  investigating  different  data  collection 
strategies  and  optimizing  these  procedures. 

2)  Improving  the  quantitative  accuracy  of  optical  images  by  exploring  the  effect  of  penalty  terms  on  the 
reconstruction  techniques.  Incorporation  of  a  priori  information  from  other  modalities  (like  MRI,  CT) 
into  the  reconstruction  procedures  and  studying  its  effect. 

3)  Exploring  effective  ways  of  displaying  and  coregistering  3D  DOT  images. 

During  last  2  years  of  funding  of  this  project,  the  main  aims  of  this  project  were  completed.  Optimization  of 
critical  computational  aspects  in  NIR  imaging  was  completed.  An  optimal  data-collection  strategy  especially  for 
the  DOT  clinical  system  at  Dartmouth  for  the  current  estimation  of  breast  tissue  optical  properties  was  also 
found.  An  effective  way  for  usage  of  a  priori  structural  of  information  from  MRI/CT  into  the  image 
reconstruction  procedure  was  developed  and  proven  that  the  quantitative  accuracy  of  DOT  images  can  be 
improved  by  at  least  a  factor  of  two  with  this  additional  information.  A  generalized  estimation  procedure  was 
developed  which  will  take  into  account  the  noise  characteristics  of  instruments  and  breast  tissue  optical 
properties  and  has  been  shown  robust  to  highly  noisy  data.  A  computationally  efficient  approach  to  dramatically 
reduce  the  size  of  the  matrix  to  be  inverted  in  cases  where  the  number  of  imaging  parameters  are  much  larger 
than  the  boundary  data  is  developed.  This  algorithm  reformulates  the  inversion  approach,  within  most  least- 
squares  approaches,  to  allow  the  inversion  to  be  based  upon  the  size  of  the  data,  rather  than  the  number  of 
imaging  field  parameters,  and  this  is  useful  for  most  3D  imaging  situations.  This  algorithm  improved  the 
computational  speed  by  at  least  a  factor  of  three.  The  final  aim  of  the  project  was  modified,  as  there  are 
commercial  platforms  available  to  coregister  3D  DOT  images  with  MRI/CT  (for  example,  OMIRAD  from 
Siemens  [4];  available  to  NIR  imaging  group  from  Siemens  as  a  part  of  network  effort  funded  by  NIH  (U54 
CA1 05480:  Network  for  Translational  Research  in  Optical  Imaging)),  to  an  important  problem  to  find  methods 
to  counteract  the  noisy  data  (which  is  more  common  in  imaging  large  breasts  (>  8  cm  in  diameter)).  This  noisy 
data  causes  large  artifacts  when  used  in  the  traditional  reconstruction  methods,  some  times  these  traditional 
methods  may  not  be  able  to  converge  to  any  meaningful  solutions  in  these  cases.  When  the  experimental  system 
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was  characterized  to  provide  noise  characteristics  of  the  data  collected,  along  with  structural  and  spectral  prior 
information,  using  these  (noise  characteristics)  in  the  generalized  reconstruction  approach  was  proven  to  give 
more  quantitative  accurate  images  (with  in  20%  error  limit)  of  spectrally-derived  tissue  functional  properties 
even  in  cases  of  highly  noisy  experimental  phantom  data.  With  out  these  noise  characteristics,  traditional 
reconstruction  algorithms  are  not  able  to  provide  any  meaningful  functional  images. 

Body 

Optimizing  the  critical  computational  aspects  of  near  infrared  tomographic  imaging 

The  image  resolution  and  contrast  in  Near- 
Infrared  (NIR)  tomographic  image 
reconstruction  are  affected  by  parameters 
such  as  the  number  of  boundary 
measurements,  the  mesh  resolution  in  the 
forward  calculation  and  the  reconstruction 
basis.  The  magnitude  of  the  total  sensitivity 
was  analyzed  to  find  the  spatial  variation 
for  a  given  problem,  and  the  field  response 
of  the  domain  becomes  more  uniform  by 
increasing  the  sensitivity  to  deeper  regions, 
while  suppressing  the  hypersensitivity  near 
the  external  boundaries.  This  is  achieved 
with  an  increase  in  the  number  of 
measurements. 

Using  singular-value  decomposition  (SVD) 
and  example  reconstructed  images, 
numbers  of  16  or  24  fibers  are  sufficient  for 
imaging  the  2D  domain.  The  number  of 
useful  measurements  actually  decreases 
exponentially  with  the  number  of 
measurements  used,  and  the  number  of 
useful  singular  values  increases  only  as  the 
logarithm  of  the  number  of  measurements. 

For  this  2D  reconstruction  problem,  given  a 
computational  limit  of  10  sec  per  iteration, 
leads  to  choice  of  forward  mesh  with  1785 
nodes  and  reconstruction  (pixel)  basis  of 
30x30  elements. 

The  use  of  three  fundamentally  different  data  collection  strategies  for  three-dimensional  (3D)  NIR  tomography 
was  compared.  Given  a  3D  NIR  imaging  problem,  using  a  single  plane  of  data  can  provide  useful  images  if  the 
anomaly  to  be  reconstructed  is  within  the  measurement  plane.  However,  if  the  location  of  the  anomaly  is  not 
known,  3D  data  collection  strategies  are  very  important.  The  recovered  quantitative  accuracy  of  the  anomaly 
region  decreases  (approximately  10%)  with  the  addition  of  out-of-plane  data  relative  to  in-plane  data.  Usage  of 
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Figure  2:  The  reconstructed  absorption  coefficient  distribution  for 
the  cylindrical  object  with  a  spherical  absorption 
inhomogeneity  (diameter  of  15mm  and  contrast  2:1  with  respect 
to  background)  located  at  x,  y  and  z  locations  (a)  (0,0,0), 

(b)  (30,0,0)  and  (c)  (30,0,10).  The  three  columns  of  images  show 
the  results  achieved  with  the  three  different  data  collection  schemes 
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single-plane  of  data  gives  slightly  better  quantitative  accuracy,  if  the  anomaly  lies  in  the  data  acquisition  plane. 
Further  the  quantitative  accuracy  of  the  reconstructed  anomaly  increased  approximately  from  15%  to  89%  as  the 
anomaly  moved  from  the  centre  to  boundary,  respectively.  The  data  supports  the  idea  that  the  use  of  in-plane 
data  in  the  3D  data  collection  strategies  may  be  sufficient  for  the  3D  NIR  imaging. 

Complete  work  along  with  the  methods  employed  and  detailed  discussion  of  the  results  given  in  the  appendix 
[5]  (Yalavarthy  et  al,  Opt.  Express  14,  p.  6113-6127,  2006). 

Usage  of  structural  a-priori  information 


NIR  tomography  combined  with  conventional 
imaging  modalities  (MRI,  CT  and  Ultra  Sound) 
has  been  a  very  active  area  of  research.  These 
hybrid  systems  show  superior  performance  in 
terms  of  qualitative  (resolution)  and  quantitative 
accuracy  compared  to  stand-alone  systems.  But 
still  there  is  lot  of  ambiguity  in  utilizing  the 
spatial  information  from  these  high  spatial 
resolution  images  into  NIR  tomography 
(coregistration).  This  work  develops  a  simple 
framework  to  incorporate  structural  a-priori 
information.  Simple  weight  matrices  that  have 
Laplacian-type  or  Helmholtz-type  structures  that 
are  derived  from  a-priori  information  have  been 
developed.  It  has  been  shown  that  utilization  of 
structural  information  using  these  weight  matrices 
will  not  bias  the  reconstruction  problem  towards  imperfect  structural  priors.  Usage  of  imperfect  a-priori 
information  in  a  parameter  reduction  (i.e.  hard-priors)  in  the  imaging  field  through  the  enforcement  of  spatially 
explicit  regions  gives  erroneous  results.  In  the  phantom  experiments,  it  is  shown  that  the  Helmholtz  type  of 
regularization  matrix  gives  the  best  estimate  of  the  scattering  parameter  and  the  Laplacian  provides  best 
estimate  for  the  absorption  parameter.  Overall,  usage  of  structural-priors  improve  the  reconstructed  image 
quantitative  accuracy  by  at  least  a  factor  of  two. 


Details  of  the  implementation  along  with  analysis 
of  results  given  in  the  appendix  [6]  (Yalavarthy  et 
al,  Opt.  Express  15,  p.  8043-8058,  2007). 

Generalized _ Least-Squares _ (GLS) 

minimization:  Two-dimensional  imaging 

(Numerical  Study) 

DOT  involves  recovery  of  the  distribution  of 
optical  parameters  by  matching  the  experimental 
data  with  modeled  data  (Levenberg-Marquardt 
(LM)  Minimization).  A  variation  of  this  approach 
by  adding  the  parameter  field  to  the  minimization 
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Figure  4:  Reconstruction  results  using  different  minimization  techniques  with  3% 
noise  in  the  data. 
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Figure  3:  Reconstruction  results  with  the  usage  of  imperfect  structural  priors  using 
different  reconstruction  techniques. 
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function  is  done  through  Tikhonov  regularization,  where  the  regularization  parameter  is  chosen  to  overcome  the 
ill-conditioning  of  the  problem.  In  this  work,  a  generalized  framework  for  DOT  was  developed  including 
variance  of  the  data  and  parameters  as  weight  matrices.  These  weight  matrices  can  also  include  the  structural 
information  obtained  by  MRI,  Ultrasound  or  X-ray  imaging.  These  weight  matrices,  include  the  system  noise 
characteristics  and  expected  size  of  optical  parameters  and  constraints  for  the  imaging  problem  and  make  the 
inversion  routine  more  robust  to  noise.  This  also  makes  the  imaging  problem  more  stable.  It  is  also  important  to 
note  that  Tikhonov  regularization  becomes  a  special  case  of  the  Generalized  Least-Squares  (GLS)  formulation. 
This  GLS  estimation  of  optical  properties  has  been  shown  to  be  very  robust  to  noise  and  proven  to  be  stable 
over  iterations. 

Complete  formulation  along  with  results  is  given  in  the  appendix  [6]  (Yalavarthy  et  al,  Med.  Phys.  34,  p.  2085- 
2098,  2007). 
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Figure  5:  Actual  and  reconstructed  (a)  pa  and  (b)  p’s  distributions  of  a  cylindrical 
target  using  experimental  multi-layered  phantom  data.  Two-dimensional  cross- 
sections  of  the  3D  volume  in  2.5  mm  increments  spanning  from  z  =  -12.5  mm  to 
z  =  12.5  mm  (from  left  to  right)  are  shown.  Actual  distributions  are  given  in  the 
first  row.  Reconstructed  distribution  using  the  Levenberg-Marquardt  (LM) 
minimization  scheme  and  GLS  minimization  scheme  are  presented  in  the  middle 
and  last  rows  respectively. _ 


A  computationally  efficient  algorithm  for 

three-dimensional  imaging 

Three-dimensional  (3D)  diffuse  optical 
tomography  (DOT)  is  known  to  be  a  non-linear, 
ill-posed  and  sometimes  under-determined 
problem,  where  regularization  is  added  to  the 
minimization  to  allow  convergence  to  a  unique 
solution.  In  this  work,  a  generalized  least- 
squares  (GLS)  minimization  method  was 
implemented,  which  employs  weight  matrices 
for  both  data-model  misfit  and  optical  properties 
to  include  their  variances  and  covariances,  using 
a  computationally  efficient  scheme.  This  allows 
inversion  of  a  matrix  that  is  of  dimension 
dictated  by  the  number  of  measurements, 
instead  of  by  the  number  of  imaging  parameters. 
This  increases  the  computation  speed  up  to  four 
times  per  iteration  in  most  of  under-determined 
3D  imaging  problems.  An  analytic  derivation, 
using  the  Sherman-Morrison-Woodbury 
identity,  is  shown  for  this  efficient  alternative 
form  and  it  is  proven  to  be  equivalent,  not  only 
analytically,  but  also  numerically.  Equivalent 
alternative  forms  for  other  minimization 
methods,  like  Levenberg-Marquardt  (LM)  and 
Tikhonov,  are  also  derived.  3D  reconstruction 
results  indicate  that  the  poor  recovery  of 
quantitatively  accurate  values  in  3D  optical 
images  can  also  be  a  characteristic  of  the 
reconstruction  algorithm,  along  with  the  target 
size.  Interestingly,  usage  of  GLS  reconstruction 
methods  reduces  error  in  the  periphery  of  the 
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image,  as  expected,  and  improves  by  20%  the  ability  to  quantify  local  interior  regions  in  terms  of  the  recovered 
optical  contrast,  as  compared  to  LM  methods.  Characterization  of  detector  (PMTs)  noise  have  enabled  the  use 
of  the  GLS  method  for  reconstructing  experimental  data  and  showed  a  promise  for  better  quantification  of  target 
in  3D  optical  imaging.  Use  of  these  new  alternative  forms  becomes  effective  when  the  ratio  of  the  number  of 
imaging  property  parameters  exceeds  the  number  of  measurements  by  a  factor  greater  than  2. 

Complete  formulation  along  with  the  analytic  derivation  is  given  in  the  appendix  [7]  (Yalavarthy  et  al,  Med.  P 
hys.  2007  (submitted)). 

Incorporation  of  noise  characteristics  into  the  near  infrared  tomographic  image  reconstruction 

algorithms 


1  1 


-3  -2.5  -2  -1.5  -i  -0.5  G 

Measured  PMT  voltage  (V.  in  volts)  in  log  scale 


Measured  PMT  voltage  [V.  in  volts)  in  log  scale 

Figure  6:  An  error  (deviation,  a)  plot  of  the  measured  voltage  and 
phase  (0)  as  a  function  of  mean  of  measured  PMT  voltage.  The  leger 
of  the  figure  represents  the  fitting  model  used.  Each  data  point 
corresponds  to  a  sample  size  of  200. 


To  incorporate  the  noise  characteristics  into  the  image 
reconstruction  algorithm,  the  first  step  was  to  characterize 
the  experimental  data  collection  system.  As 
photomultiplier  tubes  (PMT)  are  used  as  a  detectors  in  the 
data-collections  syst 

em  at  Dartmouth  for  the  data-collection  [8],  to  measure 
the  deviation  in  the  measured  signal,  a  series  of  light 
signal  measurements  were  taken  through  a  homogeneous 
intralipid  phantom  experiments  were  conducted  with 
increasing  levels  of  blood  (HbT)  concentration,  varying 
from  7.3  pM  to  36  pM,  leading  to  a  decrease  in  the 
measured  PMT  voltage.  To  achieve  this,  the  gain  of  the 
PMT  was  kept  at  0.9.  A  concise  discussion  of  the  PMT 
gain  setting  in  the  system  characterization  is  given  in  Ref. 
[8].  A  single  source  and  the  farthest  detector  (number-8) 
was  used  for  these  transmission  measurements.  For  every 


concentration  200  data  points  were  collected  to  estimate  the  deviation  in  the  measured  voltage  using  the  same 
gain  settings.  The  approach  for  the  characterization  is  similar  to  the  one  described  in  Ref.  [8],  except  the  raw 
detected  voltage  was  used  here  for  estimation  of  the  error  (or  deviation  a).  Note  also  that  two  sets  of  diameters, 
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56  mm  and  84  mm,  were  used  to  get  the  voltage  in  the  range  of  0-1  volts.  This  was  repeated  for  all  the 
wavelengths  to  ensure  uniformity  of  performance  in  the  signal,  and  to  ensure  that  the  trend  was  independent  of 
wavelength  and  gain  setting. 


Figure  6  gives  a  plot  of  error  (g(V)/V)  as  a  function  of  measured  PMT  voltage  for  785  nm  wavelength.  A 
similar  trend  was  observed  for  other  wavelengths.  This  plot  also  gives  a  deviation  in  phase  (o(0))  in  degrees  for 
the  same  voltage.  The  lowest  voltage  that  was  measured  was  0.001,  which  is  in  the  noise  floor.  The  measured 
deviations  were  1%  for  InA  (amplitude)  and  0.5°  for  0  (phase).  These  values  are  similar  to  the  ones  reported  in 
the  literature  (Mcbride  et.  al  [8]  reported  0.32%  for  the  PMT  voltage  and  0.48°  in  phase).  A  solid  line  shows 
these  average  deviation  using  1/V  fitting  (following  shot-noise  model)  in  Fig.  6.  Using  this  plot,  the  noise 
characteristics  can  be  estimated  using  the  measured  PMT  voltage.  Inclusion  of  these  noise  characteristics  into 
the  generalized  least-squares  (GLS)  scheme  combined  with  spectral  and  spatial  priors  lead  to  the  reasonable 
estimates  of  target  functional  properties  in  the  case  of  noisy  phantom  data  (example  data  set  from  a  cylindrical 
gelatin  phantom  (Fig.  7)  is  plotted  in  Fig.  8). 
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All  Det 

HbT(juM) 

sto2  (%) 

H20  (%) 
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b 
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Actual 

- 

11 

27 

- 

- 

95.0 

95.0 

- 

- 

- 

- 

LM 

NO 

10.5 

22.5 

78.1 

81.8 

86.6 

95.7 

0.71 

0.74 

0.90 

1.10 

GLS 

NO 

10.2 

22.4 

80.1 

81.9 

88.1 

92.8 

0.91 

0.89 

1.05 

1.16 

GLS 

YES 

10.4 

21.9 

81.6 

80.5 

91.6 

90.5 

0.93 

0.85 

1.08 

1.11 

Table-1:  Summary  of  reconstruction  results  of  phantom  data  (schematic  diagram  is  given  in  Fig.  7)  using  LM  and  GLS  methods.  The 
middle  two  rows  are  reconstruction  results  from  the  data  in  which  noisy  detectors  (1,2,  10,  14,  and  15)  data  points  were  removed.  The 
last  row  represents  the  case  of  GLS  scheme  using  all  the  detectors  data  using  noise  modle,  1/V2  fit  (legend  of  Fig.  6).  Notations :  FIbT: 
Total  Flemoglobin  (in  pM);  St02:  Oxygenation  Saturatin  (in  %);  H20:  Waiter  content  (in  %);  a:  Scattering  Amplitude;  b:  Scattering 
Power 


Typically  in  these  cases  of  noisy  data,  an  experienced  researcher  looks  at  the  data  and  discards  the  noisy  data 
points  from  the  reconstruction  scheme,  which  becomes  very  subjective.  The  aim  was  to  include  to  the  noisy 
characteristics  into  the  reconstruction  scheme  using  Fig.  6.  When  this  was  performed  using  all  the  data  points, 
the  estimates  using  GLS  reconstruction  scheme  was  with  in  the  20%  error  limit  of  the  expected  functional 
properties.  In  this  case  of  using  all  the  data  points,  traditional  reconstruction  algorithm  (Levenberg-Marquardt 
(LM))  did  not  give  any  meaningful  solutions.  These  results  are  summarized  in  Table- 1.  The  important 
conclusion  from  this  phantom  study  is  that  the  GLS  procedure  (developed  as  part  of  this  project)  is  able  to  give 
reasonable  estimates  of  the  functional  properties  even  when  the  data  from  the  noisy  detectors  is  used.  Where  as 
the  traditional  algorithms  are  not  able  to  converge  to  any  meaningful  solution  when  the  noisy  detectors  data  was 
used,  and  a  subjective  way  of  removing  the  noisy  detectors  data  was  able  provide  reasonable  estimates. 

All  these  results  along  with  extensive  simulation  studies  are  presented  in  the  chapter-6  of  Ref.  [9].  We  are 
panning  to  submit  a  journal  paper  with  these  results  and  analysis. 

Key  research  accomplishments 

•  Optimization  of  computational  aspects  of  DOT  image  reconstruction  (especially,  data-collection 
strategy)  {Aim-1). 

•  Effective  way  of  incorporating  structural  priors  into  NIR-DOT  image  reconstruction  procedure  (part  of 
Aim-2) 

•  Generalized  least  squares  minimization  formulation  and  extensive  testing  of  the  algorithm  in  simulations 


9 


( part  of  Aim- 2). 

•  Development  of  a  computationally  efficient  algorithm  for  under-determined  three-dimensional  DOT 
problems  (part  of  Aim-2). 

•  Inclusion  of  noise  characteristics  of  the  experimental  data  in  the  Generalized  least  squares  algorithm  to 
improve  the  estimation  capability  of  3D  functional  images  in  the  case  of  highly  noisy  data  (modified 
Aim-3). 

Reportable  outcomes 

All  the  formalities  (including  thesis  defense)  for  the  award  of  the  Ph.D.  degree  for  the  P.I.  is  completed 
(Appendix  contains  the  completion  certificate),  where  the  Ph.D.  degree  will  be  awarded  at  the  commencement 
on  June  8,  2008.  Due  to  this  training,  P.I.  is  employed  as  a  post-doctoral  research  associate  in  the  division  of 
medical  physics,  Department  of  Radiation  Oncology,  Washington  University  School  of  Medicine  in  St.  Louis. 

This  funding  has  lead  to  7  journal  publications  (including  4  first  author  publications  from  the  P.I.),  6  refereed 
international  conference  proceeding  (including  2  first  author  proceeding  from  the  P.I.),  and  5  international 
conference  presentations  (including  3  presented  by  the  P.I.).  These  are  listed  in  separate  categories  in 
chronological  order  below. 

JourimL  Publications: 

1.  P.  K.  Yalavarthy,  D.  R.  Lynch,  B.  W.  Pogue,  H.  Dehghani,  and  K.  D.  Paulsen,  “Implementation  of  a 
computationally  efficient  least-squares  algorithm  for  highly  underdetermined  three-dimensional  diffuse  optical 
tomography  problems,”  Med.  Phys.  2007  (submitted). 

2.  M.  Eames,  B.W.  Pogue,  P.  K.  Yalavarthy,  and  H.  Dehghani,  “An  efficient  Jacobian  reduction  method  for 
diffuse  optical  image  reconstruction,”  Opt.  Express  15,  15908-15919  (2007). 

3.  S.  Srinivasan,  B.  W.  Pogue,  C.  M.  Carpenter,  P.  K.  Yalavarthy,  and  K.  D.  Paulsen,  “A  boundary  element 
approach  for  image-guided  near-infrared  absorption  and  scatter  estimation,”  Med.  Phys.  34  4545 — 4557  (2007). 

4.  P.  K.  Yalavarthy,  B.  W.  Pogue,  H.  Dehghani,  C.  M.  Carpenter,  S.  Jiang,  and  K.  D.  Paulsen,  “Structural 
information  within  regularization  matrices  improves  near  infrared  diffuse  optical  tomography,”  Opt.  Express  15, 
8043-8058  (2007). 

5.  A.  L.  Darling,  P.  K.  Yalavarthy,  M.  M.  Doyley,  H.  Dehghani,  and  B.  W.  Pogue,  “Simulated  interstitial  fluid 
pressure  in  soft  tissue  as  a  result  of  externally  applied  contact  pressure,”  Phys.  Med.  Biol.  52,  4121 — 4136 
(2007). 

6.  P.  K.  Yalavarthy,  B.  W.  Pogue,  H.  Dehghani,  and  K.  D.  Paulsen,  “Weight-Matrix  Structured  Regularization 
Provides  Optimal  Generalized  Least-Squares  Estimate  in  Diffuse  Optical  Tomography,”  Med.  Phys.  34,  2085- 
2098  (2007). 

7.  P.  K.  Yalavarthy,  H.  Dehghani,  B.  W.  Pogue,  and  K.  D.  Paulsen,  “Critical  computational  aspects  of  near 
infrared  circular  tomographic  imaging:  Analysis  of  measurement  number,  mesh  resolution  and  reconstruction 
basis,”  Opt.  Express  14,  61 13-6127  (2006). 

Refereed  Conference  Proceeding: 

1.  P.  K.  Yalavarthy,  R.  Langoju,  B.  W.  Pogue,  H.  Dehghani,  A.  Patil,  and  K.  D.  Paulsen,  “Cramer-Rao 
estimation  of  error  limits  for  diffuse  optical  tomography  with  spatial  prior  information,”  Proc.  of  SPIE  6434 
(BiOS-2007  in  Photonics  West-2007,  20-25  January  2007,  San  Jose,  California),  643403:1-13  (2007). 

2.  S.  C.  Davis,  H.  Dehghani,  P.  K.  Yalavarthy,  B.  W.  Pogue,  K.  D.  Paulsen,  “Comparing  two  regularization 
techniques  for  diffuse  optical  tomography,”  Proc.  of  SPIE  6434  (BiOS-2007  in  Photonics  West-2007,  20-25 
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January  2007,  San  Jose,  California),  64340X:1-12  (2007). 

3.  H.  Dehghani,  C.  M.  Carpenter,  P.  K.  Yalavarthy,  B.  W.  Pogue,  and  J.  P.  Culver,  “Structural  a-priori 
Information  in  near  infrared  optical  tomography,”  Proc.  of  SPIE  6431  (BiOS-2007  in  Photonics  West-2007,  20- 
25  January  2007,  San  Jose,  California),  64310B:l-7  (2007). 

4.  A.  Darling,  P.  K.  Yalavarthy,  H.  Dehghani,  and  B.  W.  Pogue,  “Interstitial  fluid  pressure  due  to  externally 
applied  force  in  breast  tissue,”  Proc.  of  SPIE  6431  (BiOS-2007  in  Photonics  West-2007,  20-25  January  2007, 
San  Jose,  California),  643 10Z:  1-10  (2007). 

5.  P.  K.  Yalavarthy,  C.  Carpenter,  S.  Jiang,  H.  Dehghani,  B.  W.  Pogue,  and  K.  D.  Paulsen,  “Incorporation  of 
MR  structural  information  in  diffuse  optical  tomography  using  Helmholtz  type  regularization,”  Proc.  of  OSA 
Biomedical  Topical  Meetings,  OSA  Technical  Digest,  SH29:l-3,  Optical  Society  of  America,  Washington,  DC 
(2006). 

6.  C.  Carpenter,  B.  W.  Pogue,  P.  K.  Yalavarthy,  S.  Davis,  S.  Jiang,  H.  Dehghani,  and  K.  D.  Paulsen,  “Analysis 
of  3-dimensional  reconstruction  in  a  MR-guided  NIR  tomography  system,”  Proc.  of  OSA  Biomedical  Topical 
Meetings,  OSA  Technical  Digest,  SH33:l-3,  Optical  Society  of  America,  Washington,  DC  (2006). 

Conference  Presentations : 

1.  P.  K.  Yalavarthy,  B.  W.  Pogue,  and  H.  Dehghani,  “A  generalized  least-squares  minimization  method  for 
near  infrared  diffuse  optical  tomography,”  Department  of  Defense  Era  of  Hope  Meeting,  Baltimore,  Maryland 
25-28  June  2008  (submitted). 

2.  B.  W.  Pogue,  C.  M.  Carpenter,  P.  K.  Yalavarthy,  S.  C.  Davis,  J.  Wang,  and  K.  D.  Paulsen,  “Recovery  of 
hemoglobin  images  from  MR-guided  NIR  spectroscopy,”  SPIE  Medical  Imaging-2007,  San  Diego,  California, 
17-22  February  2007  {Oral presentation). 

3.  B.  W.  Pogue,  C.  M.  Carpenter,  P.  K.  Yalavarthy,  H.  Dehghani,  S.  Jiang,  X.  Wang,  W.  A.  Wells,  C.  A. 
Kogel,  S.  P.  Poplack,  J.  B.  Weaver,  and  K.  D.  Paulsen,  “Proposed  methods  to  improve  false  positive  and  false 
negative  rates  in  MR  breast  imaging,  through  combination  with  NIR  broadband  spectroscopy/tomography,” 
BiOS-2007  in  Photonics  West-2007,  San  Jose,  California,  20-25  January  2007  {Invited  oral  presentation). 

4.  P.  K.  Yalavarthy,  B.  W.  Pogue,  H.  Dehghani,  S.  Jiang,  and  K.  D.  Paulsen,  “Generalized  Least-Squares 
minimization  for  Magnetic  Resonance  guided  Diffuse  Optical  Tomography,”  BiOS-2007  in  Photonics  West- 
2007,  San  Jose,  California,  20-25  January  2007  {Oral  presentation). 

5.  P.  K.  Yalavarthy,  B.  W.  Pogue,  H.  Dehghani,  S.  Jiang,  and  K.  D.  Paulsen,  “Outline  of  the  Weighted  Least- 
Squares  minimization  for  Diffuse  Optical  Tomography,”  Network  for  Translational  Research  Optical  Imaging 
Network  (NTROI)  Retreat,  Hyatt  Regency  Newport  Beach,  CA,  June  22-24,  2006  {Poster presentation). 

Conclusions 


This  project  was  part  of  continuing  effort  to  develop  methods/algorithms  for  three-dimensional  alternative 
breast  imaging  modalities  at  Dartmouth.  Some  important  miles  stones  in  the  project  include  completing  the 
work  on  optimizing  the  NIR  data-collection  strategies  in  3D  (completing  Aim-1).  A  framework  to  incorporate 
the  spatial-priors  in  to  the  NIR  image  reconstruction  procedure  was  developed  and  also  proven  to  be  effective 
even  in  case  of  imperfect  spatial  priors,  which  is  part  of  Aim-2.  Moreover,  a  new  algorithm  that  takes  into 
account  noise  characteristics  of  the  instruments  was  developed  and  tested  extensively  in  the  simulation  studies. 
A  computational  efficient  algorithm  which  reduces  the  computational  cost  of  three-dimensional  imaging 
problem  by  at  least  a  factor  of  three  was  also  developed.  This  algorithm  was  compared  with  the  traditional 
reconstruction  algorithms  and  was  proven  to  be  effective  in  the  phantom  case  (completing  Aim-2).  Rather  than 
finding  an  effective  way  of  displaying  optical  images  (Aim-3),  which  could  be  achieved  using  a  commercial 
platform  (available  to  NIR  imaging  group  from  Siemens  as  a  part  of  network  effort  funded  by  NIH  (U54 
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CA1 05480:  Network  for  Translational  Research  in  Optical  Imaging)),  incorporation  of  noise  characteristics  into 
the  three-dimensional  optical  image  reconstruction  scheme  in  the  case  of  experimental  data  was  performed  and 
it  was  tested  using  phantom  data  leading  to  reasonable  estimates  (20%  of  expected  values)  of  functional 
properties  of  the  gelatin  phantom  (completing  the  modified  Aim-3). 
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Appendix 

The  Ph.D.  degree  completion  certification  of  the  P.I.  from  the  Registrar  is  attached.  Four  journal  papers  (one 
manuscript  and  three  published)  are  attached  which  were  a  direct  result  of  this  project,  where  P.I.  being  the  first 
author. 

1.  P.  K.  Yalavarthy,  H.  Dehghani,  B.  W.  Pogue,  and  K.  D.  Paulsen,  “Critical  computational  aspects  of  near 
infrared  circular  tomographic  imaging:  Analysis  of  measurement  number,  mesh  resolution  and  reconstruction 
basis,”  Opt.  Express  14,  6113-6127  (2006). 

2.  P.  K.  Yalavarthy,  B.  W.  Pogue,  H.  Dehghani,  and  K.  D.  Paulsen,  “Weight-Matrix  Structured  Regularization 
Provides  Optimal  Generalized  Least-Squares  Estimate  in  Diffuse  Optical  Tomography,”  Med.  Phys.  34,  2085- 


12 


2098  (2007). 

3.  P.  K.  Yalavarthy,  B.  W.  Pogue,  H.  Dehghani,  C.  M.  Carpenter,  S.  Jiang,  and  K.  D.  Paulsen,  “Structural 
information  within  regularization  matrices  improves  near  infrared  diffuse  optical  tomography,”  Opt.  Express  15, 
8043-8058  (2007). 

4.  P.  K.  Yalavarthy,  D.  R.  Lynch,  B.  W.  Pogue,  H.  Dehghani,  and  K.  D.  Paulsen,  “Implementation  of  a 
computationally  efficient  least-squares  algorithm  for  highly  underdetermined  three-dimensional  diffuse  optical 
tomography  problems,”  Med.  Phys.  2007  (submitted). 
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Critical  computational  aspects  of  near  infrared 
circular  tomographic  imaging:  Analysis  of 
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Abstract:  The  image  resolution  and  contrast  in  Near-Infrared  (NIR) 
tomographic  image  reconstruction  are  affected  by  parameters  such  as  the 
number  of  boundary  measurements,  the  mesh  resolution  in  the  forward 
calculation  and  the  reconstruction  basis.  Increasing  the  number  of 
measurements  tends  to  make  the  sensitivity  of  the  domain  more  uniform 
reducing  the  hypersensitivity  at  the  boundary.  Using  singular- value 
decomposition  (SVD)  and  reconstructed  images,  it  is  shown  that  the 
numbers  of  16  or  24  fibers  are  sufficient  for  imaging  the  2D  circular  domain 
for  the  case  of  1%  noise  in  the  data.  The  number  of  useful  singular  values 
increases  as  the  logarithm  of  the  number  of  measurements.  For  this  2D 
reconstruction  problem,  given  a  computational  limit  of  10  sec  per  iteration, 
leads  to  choice  of  forward  mesh  with  1785  nodes  and  reconstruction  basis  of 
30x30  elements.  In  a  three-dimensional  (3D)  NIR  imaging  problem,  using  a 
single  plane  of  data  can  provide  useful  images  if  the  anomaly  to  be 
reconstructed  is  within  the  measurement  plane.  However,  if  the  location  of 
the  anomaly  is  not  known,  3D  data  collection  strategies  are  very  important. 
Further  the  quantitative  accuracy  of  the  reconstructed  anomaly  increased 
approximately  from  15%  to  89%  as  the  anomaly  is  moved  from  the  centre 
to  boundary,  respectively.  The  data  supports  the  exclusion  of  out  of  plane 
measurements  may  be  valid  for  3D  NIR  imaging. 

©2006  Optical  Society  of  America 

OCIS  Codes:  (170.0170)  Medical  optics  and  biotechnology;  (100.3190)  Inverse  problems; 
(170.3660)  Light  propagation  in  tissues;  (170.4580)  Optical  diagnostics  for  medicine; 
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1.  Introduction 

In  the  recent  years,  there  has  been  a  heightened  interest  in  near-infra-red  (NIR)  optical 
tomography,  for  applications  such  as  diagnostic  breast  cancer  imaging  [1-3]  and  for  brain 
function  assay  [1,  4,  5].  In  NIR  tomography,  the  aim  is  to  reconstruct  interior  optical 
properties  of  the  tissue  under  investigation  from  a  finite,  yet  incomplete  set  of  transmission 
measurements  taken  at  the  tissue  external  boundaries.  The  reconstructed  optical  properties  can 
give  clinically  useful  information  regarding  tissue  physiology  and  state,  such  as  chromophore 
concentration  and  oxygen  saturation.  Typically,  the  optical  source  light  used  for  excitation  in 
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NIR  studies  is  delivered  through  optical  fibers  and  the  transmitted  light  is  also  collected 
through  the  same  or  additional  fibers  which  are  in  contact  with  the  external  surface  of  the 
tissue.  Using  these  measurements,  distributions  of  wavelength  dependent  absorption  and/or 
scattering  coefficients  of  the  tissue  are  reconstructed  using  a  model-based  iterative  algorithm. 
NIR  studies  have  the  advantage  of  being  non-invasive,  non-hazardous  and  can  therefore  be 
applied  repeatedly  to  investigate  functional  changes  in  tissue  over  a  prolonged  time. 

The  dominance  of  light  scattering  in  tissue  at  NIR  wavelengths  makes  optical  tomography 
inherently  more  difficult  in  the  sense  that  light  becomes  diffuse  within  millimeters  of  travel, 
reducing  the  resolution  of  the  reconstructed  images.  The  image  reconstruction  procedure  (i.e. 
the  inverse  problem)  is  non-linear,  ill-posed  and  ill-conditioned  [6]  and  to  improve  image 
reconstruction,  the  number  of  measurements  are  generally  increased,  to  increase  the  amount 
of  independent  information.  However  due  to  experimental  set-up  constraints,  such  as  the  light 
collection  strategy,  source  and  detector  fiber  size  and  the  imaging  domain  geometry,  the  total 
number  of  boundary  measurements  that  can  be  taken  from  is  often  quite  limited.  In  addition, 
there  are  constraints  on  the  data  acquisition  and  computation  time  that  need  to  be  considered 
for  the  specific  application  in  which  NIR  light  is  used. 

There  have  been  some  limited  studies  [7-11]  on  optimization  of  the  fiber  positions  and 
measurements  to  get  the  best  possible  image  resolution  and  contrast  in  NIR  tomography. 
More  specifically,  Culver  et.  al  [11]  have  showed  that  singular  value  decomposition  (SVD) 
analysis  of  the  weight  matrix  (also  known  as  the  Jacobian  or  sensitivity  matrix)  can  be  used  to 
optimize  detector  placement  in  the  reflectance  and  direct  transmittance  geometries  of  a 
homogeneous  slab  medium,  and  indicated  that  this  could  be  extended  to  arbitrary  geometries 
with  heterogeneous  tissue  volumes.  However,  there  remain  many  unknowns  regarding  the 
appropriate  number  of  measurements  required  to  get  a  sufficiently  good  image  given  the 
practical  constraints  of  measurement  number  and  image  recovery  algorithm,  which  is  the 
subject  of  this  paper.  Furthermore,  few  studies  have  specifically  investigated  the  effect  of 
mesh  resolution  in  both  the  forward  and  inverse  calculations  and  very  little  is  known  about  the 
quantitative  increase  in  accuracy  which  is  a  direct  result  of  mesh  resolution  and  appropriate 
reconstruction  bases.  This  work  is  an  attempt  to  answer  questions  regarding  the  limited 
increase  in  number  of  measurements,  more  specifically  benefits  from  the  increased  amount  of 
information  as  well  as  investigating  aspects  that  will  have  effects  on  image  reconstruction 
procedure  and  resolution  as  well  as  the  contrast  of  the  reconstructed  image. 

In  the  present  work,  both  a  two  dimensional  (2D)  circular  domain  and  a  three  dimensional 
(3D)  cylindrical  geometry  are  investigated  since  most  investigations  to  date  have  used  either 
of  these  geometries  for  system  and  algorithm  evaluation.  Initially  the  effect  of  mesh  resolution 
is  investigated  in  the  forward  problem  by  comparing  the  Jacobian  cross-section  for  various 
resolution  2D  meshes  to  show  improvements  in  numerical  accuracy.  Next  the  effect  of 
increasing  the  number  of  measurements  upon  the  resulting  reconstructed  image  using 
singular-value  analysis  is  investigated.  Results  regarding  the  optimized  reconstruction  basis 
are  presented  for  the  given  2D  model,  and  the  impact  in  the  Root  Mean  Square  (RMS)  error  of 
increased  spatial  sensitivity  is  presented  as  a  function  of  increasing  number  of  measurements. 
Finally  a  case-to-case  analysis  is  shown  by  increasing  the  number  of  measurements  in  image 
reconstruction  procedure  and  comparing  the  underlying  image  errors  within  the  reconstructed 
images. 

Since  3D  problems  have  more  degrees  of  freedom  (unknown  parameters),  they  are  highly 
ill-determined  as  compared  to  the  2D  problem.  But  NIR  optical  tomography  utilizes  the  data 
from  the  3D  tissue  volumes  and  therefore  should  be  treated  as  a  3D  imaging  problem.  Since 
light  propagation  in  tissue  is  physically  spread  in  all  directions,  3D  models  are  known  to  be  an 
accurate  prediction  of  the  light  fluence,  whereas  2D  models  are  simple  yet  inaccurate  at 
predicting  the  interior  fluence  distributions  [4,  12-17].  In  order  to  further  advance  NIR  optical 
tomography  into  a  suitable  and  accurate  clinical  imaging  modality,  it  is  important  to  develop 
fully  3D  imaging  tools,  yet,  the  major  challenge  in  this  task  is  to  determine  how  to  acquire 
large  data  sets  which  overcome  the  inherent  limitation  of  the  3D  problem  being  ill-determined 
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[18].  That  is,  to  improve  image  reconstruction  quality  in  3D,  the  number  of  measurements  can 
be  increased  as  mentioned  in  2D  case,  even  here  these  measurements  are  quite  limited. 

For  the  chosen  3D  cylindrical  geometry,  for  example,  acquiring  experimental  data  from 
three  different  planes  of  fiber  setup  improves  the  reconstructed  image  of  the  entire  domain  as 
compared  to  one  single  plane  of  data,  as  there  are  greater  numbers  of  measurements  providing 
a  larger  set  of  sampling  of  the  entire  volume  of  interest.  There  are  many  strategies  to  increase 
measurement  number  and  it  is  not  clear  which  present  the  best  improvement  in  the  final 
image.  Specifically,  this  work  examines  effects  of  different  measurement  strategies  for  3D 
NIR  tomography  by  presenting  and  quantifying  the  underlying  effects  of  using  a  single  plane 
of  tomographic  data  as  compared  to  three  planes  of  tomographic  data.  Within  the  latter  case, 
this  work  also  presents,  quantifies  and  discusses  the  benefits,  limits  and  losses  due  to  the 
measurement  of  in-plane  data  as  compared  to  out-of  plane  data  and  will  compare  and  contrast 
these  data  collection  geometries  from  the  prospective  of  gain  and  loss  in  the  reconstructed 
image  quality  and  respective  computation  time. 

2.  Methods 


Conventional  numerical  methods  for  the  forward  calculations  in  NIR  imaging  use  the  Finite 
Element  Method  (FEM),  which  is  considered  as  a  flexible  and  accurate  approach  to  modeling 
heterogeneous  domains  with  arbitrary  boundaries.  Light  transport  in  scattering  tissue  can  be 
accurately  described  by  the  Diffusion  Approximation  (DA)  to  the  Radiative  Transfer  Equation 
(RTE)  [19]: 


-  V.Ar(r)V<D(r,&>)  + 


Ma(r)  + 


O (r,co)  -  q0(r,co) 


(1) 


where  0(r,tf/)is  the  photon  density  at  position  r  and  modulation  frequency  CO  (100  MHz  in 

this  work),  and  k  =  l/[3(qa  +  |a/)],  the  diffusion  coefficient,  where  qa  and  |a/  are  the 
probabilities  per  unit  length  of  absorption  and  transport  scattering,  respectively,  and 

q0  (r,  CO)  is  an  isotropic  source  term.  The  Robin  (Type  III)  boundary  condition  is  used  which 

best  describes  the  light  interaction  from  a  scattering  medium  to  the  external  air  boundary  [20]. 
The  calculated  boundary  data  values  with  a  frequency  domain  system  are  the  amplitude  and 
phase  of  the  signal,  from  which  the  diffusion  and  absorption  coefficients  can  be 
simultaneously  reconstructed. 

For  the  inverse  problem,  a  small  change  in  boundary  data  is  related  to  a  small  change  in 
optical  properties  through  the  Jacobian  matrix  of  values.  The  Jacobian  matrix  for 
reconstructing  both  the  unknowns  using  two  different  data-types  is  calculated  using  the 
Adjoint-method  [21],  and  has  dimensions  of  (2xSxD)  by  (2xN),  where  S  and  D  are  the 
number  of  sources  and  detectors  corresponding  to  each  source  respectively.  N  represents  the 
number  of  nodes  in  the  mesh  used  in  the  forward  calculation.  Here  the  Jacobian  maps  the 
changes  in  log  amplitude  and  phase  (2xSxD)  to  both  absorption  and  diffusion  changes  at  each 
node  of  the  FEM  model  (2xN).  The  Jacobian  which  maps  the  change  in  detected  signal  to 
image  space  has  four  parts: 
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In  all  our  analysis,  only  the  J2  section  is  considered  (dimension  of  (SxD)  by  N),  which  maps  a 
small  change  in  the  absorption  coefficient  to  a  small  change  in  measured  log  intensity  of  the 
signal.  Since  all  kernels  of  the  complete  Jacobian  show  similar  results,  the  discussion  is 
limited  to  the  results  of  /2,  and  shall  henceforth  be  referred  to  as  J. 

In  the  reconstruction  procedure  presented,  a  modified  Levenberg-Marquadt  algorithm  is 
used  for  calculating  the  estimates  of  qa,  which  is  an  iterative  procedure  [10]  solving: 
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[A|i.a]  =  [JTJ  +  XI]-'.  JTb 


(3) 


Here  [A|Ha]  is  an  update  vector  for  the  absorption  coefficient,  I  is  the  identity  matrix  and  X  is  a 
regularization  parameter.  Also,  b  =  [y  -  a)],  where  y  is  the  measured  (or  simulated) 

heterogeneous  boundary  data  and  jF(pa)  is  the  forward  data  for  the  current  estimate  of  pa.  In 
all  of  the  presented  work  using  simulated  data,  1%  noise  was  added  to  the  amplitude,  which  is 
a  typical  noise  observed  in  experimental  data  [2]. 

For  the  2D  analysis  a  circular  model  with  a  diameter  of  86  mm  centered  at  (0,  0)  and  with 
homogeneous  optical  properties  of  pa  =  0.01  mm'1  and  p/  =1.0  mm'1  is  considered.  The  light 
collection/delivery  fibers  are  arranged  in  a  circular  equally  spaced  fashion,  where  one  fiber  is 
used  as  the  source  while  all  other  fibers  are  used  as  detectors,  to  give  T’  number  of 
measurements  (where  P=  M(M-l),  where  M  is  number  of  fibers).  The  source  is  a  Gaussian 
source  of  Full  Width  Half  Maximum  (FWHM)  of  3mm,  and  it  is  placed  one  transport 
scattering  length  within  the  external  boundary. 


Fig.  1 .  Schematic  diagram  of  data  collection  geometry  used  for  the  3D  cylindrical  model. 

For  the  3D  analysis,  a  cylindrical  medium  with  a  diameter  of  86  mm  having  height  of  100 
mm  centered  at  (0,  0,  0),  with  homogeneous  optical  properties  of  pa  =  0.01  mm'1  and  p/  =  1 
mm'1  is  used  (Fig.  1).  The  light  collection/delivery  fibers  are  arranged  in  a  circular  and 
equally  spaced  fashion  and  are  in  either  a  single  plane  of  16  fibers  or  3  planes  of  16  fibers  per 
plane,  totaling  48  fibers.  Specifically  three  different  strategies  for  data  collection  are 
considered: 

(a) .  Single  layer  data:  The  16  fibers  are  arranged  in  a  circular  and  equally  spaced  fashion  in  a 
single  Layer-I  (Fig.  1),  where  one  fiber  is  used  at  a  time  as  the  source  while  all  other  fibers  are 
used  as  detectors,  to  give  240  (16x15)  amplitude  measurements. 

(b) .  Three  layers  of  in-plane  data:  The  48  fibers  are  arranged  in  a  circular  equally  spaced 
fashion  in  all  three  layers  (Layers-I,  II  &  III  in  Fig.  1),  giving  16  fibers  per  plane,  where  one 
fiber  is  used  at  a  time  as  the  source  while  only  those  fibers  in  the  same  “source  fiber  layer”  are 
used  as  detectors,  to  give  720  (3x16x15)  amplitude  measurements. 

(c) .  Three  layers  of  out-of-plane:  Same  as  above,  except  when  one  fiber  is  used  at  a  time  as 
the  source,  all  other  fibers  in  all  three  planes  are  used  as  detectors.  This  leads  to  2256  (48x47) 
amplitude  measurements. 

For  the  image  reconstruction  process,  an  iterative  update  to  the  Jacobian  matrix  was 
computed,  after  each  successive  image  estimation.  At  each  iteration,  the  objective  function 
was  evaluated  to  estimate  the  projection  error.  The  reconstruction  procedure  was  then  stopped 
when  the  projection  error  decreased  by  less  than  3%. 
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Fig.  2.  The  sensitivity  (Jacobian)  contour  plot  of  log  amplitude  and  ga  for  a  source  (S)  and 
detector  (D),  which  are  diagonally  opposite  to  each  other  as  shown,  calculated  on  a  circular 
mesh  of  9664  nodes. 

2.1.  2D  Mesh  Resolution 

In  FEM  the  domain  is  divided  into  finite  discretized  sub-domains  wherein  the  numerical 
accuracy  and  stability  depends  highly  on  this  discretization  (mesh  resolution).  Since  the 
Jacobian  represents  the  sensitivity  of  the  detected  signal  to  a  small  change  in  optical 
properties,  the  numerical  accuracy  of  this  value  is  crucial  component  of  the  image 
reconstruction  problem,  to  study  the  effect  of  mesh  resolution  in  2D  case,  we  choose  different 
resolution  meshes  (with  number  nodes  ranging  from  150  to  4617  nodes)  along  with  a  high- 
resolution  mesh  of  9664  nodes  for  calculation  of  Jacobian.  The  Jacobian  with  a  diagonally 
opposite  source  and  detector  is  used,  as  shown  in  Fig.  2,  from  which  the  RMS  error  is 
calculated  for  each  mesh  with  respect  to  the  high-resolution  mesh.  The  RMS  error  is 
calculated  by  interpolating  the  Jacobian  of  each  mesh  unto  a  uniformly  distributed  grid, 
allowing  direct  comparison  of  each  result.  Since  the  Jacobian  represents  the  sensitivity  of  the 
detected  signal  to  a  small  change  in  optical  properties,  the  numerical  accuracy  of  this  value  is 
a  crucial  component  of  the  image  reconstruction  problem.  Here  the  highest  resolution  mesh 
provides  the  most  accurate  and  numerically  stable  solution,  therefore  the  calculated  RMS 
error  indicates  the  numerical  accuracy  of  each  lower  resolution  mesh.  The  computation  time 
taken  for  calculation  of  Jacobian  and  forward  data  is  also  noted  as  a  function  of  mesh 
resolution.  All  the  computations  were  carried  out  on  Pentium-IV  2.5  GHz  processor  with  2 
GB  of  RAM. 

2.2.  Singular-Value  (SV)  analysis 

Singular- Value  (SV)  analysis  for  the  Jacobian  matrix  is  explained  in  detail  elsewhere  [10]. 
Using  SV-analysis,  the  Jacobian  is  decomposed  into: 

J  =  USVT  (4) 

where,  U  &  V  are  orthonormal  matrices  containing  the  eigenvectors  of  J  and  S  is  a  diagonal 
matrix  containing  the  singular  values  of  J.  Vectors  of  U  and  V  correspond  to  the  modes  in  the 
detection  space  and  image  space,  respectively,  while  the  magnitude  of  the  singular  values  in  S 
represents  the  importance  of  the  corresponding  eigenvectors  in  U  and  V.  More  nonzero 
singular  values  indicating  more  modes  are  effective  in  between  the  two  spaces,  which  bring 
more  detail  and  improve  the  resolution  in  the  space.  There  are  normally  P  nonzero  singular 
values  in  the  diagonal  matrix  and  these  values  are  sorted  in  decreasing  order.  Typically  only 
those  singular  values  above  the  noise  level  (in  this  study,  1  %  noise  in  amplitude)  are  used,  as 
they  contain  the  only  useful  information  in  the  matrix.  Thus,  it  is  possible  to  determine 
whether  increasing  the  number  of  measurements  gives  rise  to  an  increase  in  the  number  of 
useful  singular  values,  which  indicates  improvement  in  the  recovered  images. 
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In  2D,  this  analysis  was  applied  to  two  separate  cases:  (1)  a  homogeneous  case  with  optical 
properties  as  given  before,  and  (2)  a  heterogeneous  case  which  mimics  breast  optical 
properties  [22],  with  properties  of  fibro-glandular  layer  being  pa  =  0.003  mm'1  and  p/  =  0.95 
mm'1  and  having  diameter  of  66  mm  and  fatty  layer  surrounding  it  having  pa  =  0.006  mm'1 
and  p/  =  1.1  mm'1  with  a  thickness  of  20mm.  The  number  of  useful  singular  values  above  the 
noise  level  were  calculated  as  the  number  of  measurements  was  increased.  The  mesh  that  was 
found  to  have  an  optimum  resolution  from  the  previous  analysis  of  the  Jacobian  (Sec.  2.1)  was 
used  for  these  analysis.  For  both  these  cases,  the  percentage  of  useful  measurements  with 
respect  to  total  number  of  measurements  was  calculated  as: 


Useful  measurements  (in  %)  = 


Useful  number  of  singular  values 
Total  number  of  singular  values 


xlOO 


(5) 


Additionally,  the  effect  of  mesh  resolution  was  studied  for  its  impact  on  the  number  of 
independent  boundary  data  points  with  an  increase  in  number  of  measurements  by  calculating 
the  rank  of  the  Jacobian,  which  is  defined  as  the  maximum  number  of  linearly  independent 
rows/columns  of  a  given  matrix.  As  each  row  of  the  Jacobian  indicates  each  measurement,  the 
rank  of  the  Jacobian  indicates  the  total  number  of  independent  measurements. 

Image  reconstruction  consists  of  two  separate,  yet  equally  important  parts;  the  forward 
model  and  the  inverse  model.  For  the  forward  model,  the  mesh  used  in  FEM  needs  to  be  such 
that  to  ensure  numerical  accuracy,  as  already  discussed.  For  the  inverse  problem,  however,  the 
goal  is  to  reduce  the  number  of  unknowns  for  the  iterative  update  by  the  use  of  a 
reconstruction  basis  [23].  Therefore  it  is  important  to  investigate  the  effects  of  various 
reconstruction  basis  degrees  of  freedom  on  the  reconstruction.  Various  reconstruction  basis 
can  be  used,  such  as  second  mesh  basis  [24],  pixel  basis  [23]  or  adaptive  [25,  26]  .  With  this 
goal,  a  reconstruction  basis  was  optimized  for  the  given  2D  problem  by  looking  at  the  number 
of  useful  singular  values  for  various  pixel  (reconstruction)  basis.  A  linear  pixel  basis  of  having 
100  (10  by  10)  elements  to  1600  (40  by  40)  elements  was  used  and  the  Jacobian  was  mapped 
to  this  basis  for  the  analysis. 


Table  1.  The  RMS  error  (with  respect  to  the  fine  mesh  of  9664  nodes)  in  the  Jacobian  cross-section  from  center  to 
boundary,  (indicated  by  dashed  line  in  Fig.  2)  at  y  =  0  mm.  This  is  tabulated  as  a  function  of  mesh  resolution,  or 
number  of  nodes  in  the  mesh.  Last  two  rows  show  the  computation  time  taken  for  calculation  the  Jacobian  and 
Forward  data  for  16  source- detector  pairs  (240  measurements).  For  the  fine  mesh  of  9664  nodes  the  computation  time 
for  Jacobian  and  Forward  data  is  98.1  sec  and  28  sec  respectively. 


Nodes 

150 

425 

1360 

1785 

2683 

3047 

3569 

4617 

RMS  error 

60.56 

27.84 

5.06 

4.84 

2.57 

2.15 

1.85 

1.07 

Jacobian 
Computation 
Time  (in  Sec.) 

1.1 

2.5 

7.8 

10.1 

15.2 

17.8 

20.8 

38.1 

Forward  data 
Computation 
Time  (in  Sec.) 

0.1 

0.3 

0.9 

1.2 

1.9 

2.2 

2.6 

9.8 

2.3.  Reconstruction  examples 

In  order  to  understand  the  effect  of  increasing  the  number  of  measurements  on  total  sensitivity 
for  a  given  2D  model  the  magnitude  of  the  Jacobian  was  examined  as  a  function  of  number  of 
measurements.  To  achieve  this,  the  horizontal  cross-section  of  the  whole  Jacobian  was 
plotted,  which  was  summed  up  over  all  measurements,  from  center  to  boundary,  and 
examined  as  the  number  of  measurements  increased.  Since  the  Jacobian  provides  relative 
sensitivity,  a  cross-section  plot  was  normalized  in  each  case  with  respect  to  its  magnitude  at 
the  center  of  the  model  and  calculated  as  a  function  of  number  of  measurements  (56  to  4032). 
For  the  3D  model,  the  cross-section  of  the  total  Jacobian  was  normalized  with  respect  to  its 
magnitude  at  the  center  of  the  model  (as  in  the  2D  case),  for  each  case  of  the  three  different 
data  collection  strategies.  Finally,  for  the  2D  model,  only  the  absorption  coefficient  was 
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reconstructed  with  an  increasing  number  of  measurements  of  an  object  with  absorption 
inhomogeneity  at  various  positions  of  domain  using  log  of  amplitude  data.  A  circular 
absorption  anomaly  of  diameter  of  10  mm  was  used  having  a  contrast  of  2:1  compared  to  its 
background.  We  used  the  optimal  forward  mesh  along  with  optimal  reconstruction  basis  for 
the  reconstruction  procedure.  A  total  of  2  positions  of  absorption  inhomogeneity  were 
considered  with  it  center  at  (x,y)  of  (0,  0),  and  (30,  0)  for  various  number  of  measurements 
starting  from  56  to  4032. 


Fig.  3.  Singular  value  analysis  of  homogeneous  and  heterogeneous  2D  circular  models,  (a). 
Plot  of  the  useful  singular  values  versus  number  of  measurements,  (b).  Plot  of  percentage  of 
useful  measurements  versus  the  total  number  of  measurements,  (c).  Plot  of  the  Rank  versus 
number  of  measurements  is  shown  for  a  range  of  mesh  nodes,  (d).  Plot  of  the  number  of  useful 
singular  values  versus  number  of  measurements  is  shown,  for  various  reconstruction  bases. 


For  the  3D  case,  a  spherical  absorption  anomaly  of  diameter  of  15  mm  was  assumed 
having  a  contrast  of  2:1  compared  to  its  background.  A  total  of  3  positions  of  absorption 
inhomogeneity  were  considered  with  its  center  at  x,  y  and  z  of  (0,0,0),  (30,0,0)  and  (30,0,10). 
The  anomalies  were  reconstructed  using  the  noise  added  data  (1%  in  amplitude)  simulated 
from  the  three  different  fiber  location  strategies.  Full  Width  at  Half-Maximum  (FWHM)  was 
measured  for  each  of  the  peaks  in  the  X-Y  and  Z-Y  planes  as  well  as  the  total  computation 
time  for  reconstruction  process. 


Table  2.  The  number  of  useful  measurements  above  the  1%  expected  noise  level,  is  shown  for  the  2D  circular  and  3D 
cylindrical  models,  having  16  source  and  detector  fibers  with  one  or  three  planes  of  data  collection.  The  two  upper 
rows  have  only  1  plane  of  collection,  whereas  the  2nd  last  row  has  3  planes  of  collection  but  not  between  the  planes, 
and  the  last  row  has  3  planes  of  data  collection  with  complete  out  of  plane  measurements. 


Number  of 
Unknowns 

Number  of 
Measureme 
nts 

Number  of 
Useful  Singular 
values 

Useful 

measurements 

(%) 

Magnitude  of  largest 
singular  value 

2D 

1785 

240 

91 

37.92 

796.4 

3D  1  layer 

20163 

240 

107 

44.58 

117.1 

3D  3layer  in¬ 
plane 

20163 

720 

269 

37.36 

164 

3D  3layer  out- 
of-plane 

20163 

2256 

328 

14.54 

304.6 

3  Results 

Figure  2  shows  a  sensitivity  plot  of  log  amplitude  and  the  absorption  coefficient  using  a  2D 
mesh  with  9664  nodes  for  a  source  and  detector  which  are  diagonally  opposite  to  each  other. 
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Table- 1  shows  the  RMS  error  with  respect  to  the  high  resolution  mesh  in  the  horizontal  cross- 
section  (as  indicated  by  the  dotted  line  in  Fig.  2)  using  the  method  described  earlier.  The  RMS 
error  calculated  here  was  also  calculated  along  different  cross-sections  of  the  model  and  a 
similar  trend  was  seen.  The  mesh  with  1785  nodes  was  found  to  have  an  RMS  error  of  less 
then  5%  as  compared  to  the  finest  mesh. 


Fig.  4.  Comparison  of  Jacobian  cross-section  with  respect  to  measurement  number,  (a).  The 
horizontal  cross-sectional  plot  of  the  sum  of  2D  circular  Jacobian  matrix  values,  from  center  to 
the  boundary  at  y  =  0  mm.  (b)  The  normalized  sum  of  2D  circular  Jacobian  matrix  values,  with 
respect  to  the  value  at  the  center  (at  x  =  0  mm,  y  =  0  mm).  The  legend  gives  number  of 
measurements  associated  with  each  plot. 

The  2D  mesh  with  1785  nodes  was  used  for  the  calculation  of  the  Jacobian  and  the 
expected  noise  level  in  the  amplitude  measurements  was  assumed  to  be  1%.  For  both  the 
heterogeneous  and  homogeneous  2D  cases,  the  number  of  useful  singular  values  above  the 
noise  level  were  calculated,  and  the  results  are  shown  in  Fig.  3(a).  Figure  3(b)  is  a  bar  chart 
showing  useful  measurements  in  percentage  [given  by  Eq.  (5)]  for  each  set  of  measurements. 
Figure  3(c)  is  a  plot  of  the  rank  of  the  Jacobian  versus  the  total  number  of  measurements  for 
meshes  having  different  resolution  starting  from  150  to  3569  nodes  versus  number  of 
measurements.  The  Jacobian  calculated  is  also  mapped  onto  a  reconstruction  (pixel)  basis 
ranging  from  10  x  10  to  40  x  40.  The  number  of  useful  singular  values  as  function  of  pixel 
basis  elements,  for  each  set  of  measurements,  are  plotted  in  Fig.  3(d).  Finally,  for  the  2D  case, 
Fig.  4  shows  the  total  sensitivity  distribution  at  the  mid-axis  cross-section,  as  a  function  of  the 
number  of  measurements.  Table  2  shows  the  number  of  useful  singular  values  of  the  3D 
model  Jacobian  which  are  above  the  noise  level  (1%)  for  the  three  different  strategies,  and 
indicates  the  effective  number  of  measurements  which  will  be  contributing  to  the 
reconstructed  image  space  and  quality.  The  number  of  useful  singular  values  is  higher  for  the 
three  layer  out-of-plane  strategy.  The  useful  percentage  of  measurements  is  higher  for  the  3D 
single  plane  of  data,  whereas  the  condition  number  is  very  high  for  the  3D  three-layer  out  of 
plane  case.  Similar  data  is  also  included  using  the  2D  circular  geometry  for  comparison 
purposes,  with  240  measurements  and  the  same  corresponding  optical  properties  as  the  3D 
model. 

The  plots  of  the  3D  Jacobian  magnitude  as  normalized  to  the  value  at  the  center  of  the 
model  are  shown  in  Fig.  5.  These  plots  shows  that,  all  the  three  strategies  of  data  collection  in 
3D  are  hypersensitive  (in  X  &  Y  direction)  at  the  boundary.  Moreover  this  is  pronounced  for 
the  3D  single-plane  case.  In  the  Z-direction  (not  shown)  it  was  found  that,  as  expected  that, 
the  sensitivity  decreases  as  the  position  moves  from  centre  to  boundary  for  all  the  three  cases. 
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Z=  0  mm 


Z=  10  mm 

Fig.  5.  The  normalized  cross-section  in  the  X-' 
dotted  line  in  Fig.  2,  from  x=  0  mm  to  x  =  43  m 
with  respect  to  the  sensitivity  at  the  origin,  (i.e. 


Z=  5  mm 


Z=  15  mm 

plane,  showing  the  total  sensitivity  across  the 
(center  to  boundary)  at  Y  =  0  mm  normalized 
X  =  0,  Y  =  0  &  Z  =  0  mm). 


Original  (la  56  240  552  992 


1560  2256  3080  4032 


0.01  0.012  0.014  0.016  0.018  0.02 


Fig.  6.  The  reconstruction  of  the  ga  distribution,  using  noisy  simulated  data  of  log  amplitude, 
for  a  circular  object  with  an  absorbing  inhomogeneity  at  the  center.  Different  numbers  of 
measurements  were  used  as  denoted  above  each  image,  ranging  from  56  up  to  4032  data  points. 

The  forward  mesh  was  1785  nodes  and  the  pixel  basis  consisted  of  30x30  elements.  The 
original  ga  distribution  is  shown  as  the  first  image. 

The  2D  reconstruction  of  a  circular  object  with  a  centralized  absorption  anomaly  of 
diameter  of  10mm  using  different  number  of  measurements,  along  with  original  pa 
distribution,  is  shown  in  Fig.  6.  The  contrast  of  the  inhomogeneity  to  background  is  2:1  and 
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for  these  reconstructions  a  pixel  basis  of  30  x  30  elements  is  used,  with  a  forward  mesh 
consisting  of  1785  nodes.  Figure  7  shows  the  plot  of  logarithm  of  rms  error  in  the  horizontal 
cross-section  (as  sown  by  dotted  line  in  Fig.  6)  as  a  function  of  measurement  number.  The 
legend  of  the  figure  gives  the  position  of  the  inhomogeneity  (diameter  of  10mm). 

Table  3  summarizes  the  results  of  the  3D  reconstruction.  Figure  8(a)  shows  the 
reconstructed  absorption  coefficient  distributions  for  a  spherical  absorption  inhomogeneity 
(diameter  of  15  mm)  located  at  (0,  0,  0)  with  a  contrast  of  2:1  to  background,  using  the  data 
collected  from  the  three  strategies.  Figure  8(b)  shows  the  results  of  the  same  effort  with  a 
spherical  inhomogeneity  located  near  to  the  boundary  (30,  0,  0).  The  results  show  that  the 
quantitative  values  of  the  anomaly  increases  as  the  anomaly  is  moved  from  centre  to  boundary 
in  X  &  Y  direction.  The  anomaly  for  this  location  is  reconstructed  with  89%  quantitative 
accuracy  compared  to  the  15%  accuracy  for  central  location.  Finally  the  reconstructed 
absorption  coefficient  distribution  for  a  spherical  absorption  inhomogeneity  (diameter  -  15 
mm),  which  is  centered  at  (30,  0,  10)  are  shown  in  Fig.  8(c)  and  it  can  be  seen  that  single  layer 
case  reconstructed  the  anomaly  in  the  wrong  location.  Here,  both  the  in-plane  and  out-of- 
plane  strategies  are  able  to  give  up  to  84%  quantitative  accuracy  (Table  3). 


Fig.  7.  A  plot  of  logarithm  of  rms  error  in  the  horizontal  cross-section  of  |la  at  y  =  0  (as  shown 

in  original  |la  distribution  of  Fig.  6)  versus  number  of  measurements  for  various  positions  of  an 
absorption  inhomogeneity.  These  calculations  used  1785  nodes  in  the  mesh  of  the  forward 
problem  and  a  pixel  basis  of  30x30  elements  in  the  reconstruction. 

4  Discussion 

The  decrease  in  the  RMS  error  for  the  horizontal  cross-section  of  the  2D  Jacobian  for  a  given 
source-detector  (diagonally  opposite  each  other)  for  a  mesh  greater  than  1500  nodes  as 
compared  to  9664  nodes  (Table- 1)  is  below  5%.  It  should  be  noted  that  the  other  kernels  of 
the  Jacobian,  for  example  J3  (d<9),  showed  better  accuracy  (2%)  when  the  mesh  had  1785 

die 

nodes  or  greater.  As  with  many  iterative  reconstruction  problems,  optical  tomography  requires 
repeated  forward  calculations  and  re-computation  of  the  Jacobian,  thereby  increasing  mesh 
resolution  which  further  implies  increase  in  computational  time,  which  is  clearly  evident  from 
last  two  rows  of  Table  1.  A  computation  limit  of  10  seconds  per  iteration,  lead  to  a  choice  of 
mesh  resolution  with  1785  nodes  for  the  forward  problem  in  two-dimensional  case,  and 
extending  this  same  level  of  resolution  to  3D  would  require  nearly  80,000  nodes,  which  is 
near  the  limit  of  what  can  be  done  computationally.  Thus  much  of  the  2D  study  presented 
here  was  run  at  the  level  of  1785  nodes.  Since  the  computation  of  the  Jacobian  using  the  FEM 
relies  on  the  discretization  of  the  domain  and  the  accuracy  of  the  numerical  model  depends  on 
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Table  3.  The  computation  time  and  accuracy  of  the  3D  reconstruction  is  shown  for  the  three 
different  data  collection  strategies,  along  with  three  different  locations  of  the  anomaly  for  each. 


Strategy 

Position  of 
anomaly 
(original) 

Iterati 

ons 

Total 

Computation 
time  (s) 

Quantitative 
accuracy  (%)  of 
the  reconstructed 
anomaly 

FWHM 
along  X- 
axis  (mm) 

FWHM 

along 

Z-axis 

(mm) 

3D:  1  layer 

(0,0,0) 

11 

3179 

15% 

16.1 

25.2 

(30,0,0) 

14 

4046 

89% 

17.2 

23.3 

(30,0,10) 

10 

2890 

- 

- 

- 

3D  3layer  in- 
plane 

(0,0,0) 

14 

8022 

14% 

16.5 

25.3 

(30,0,0) 

14 

8022 

80% 

13.1 

18.7 

(30,0,10) 

12 

6876 

110% 

11.2 

18.6 

3D  3layer 
out-of-plane 

(0,0,0) 

6 

10926 

11% 

23.7 

24.1 

(30,0,0) 

9 

16389 

78% 

13.6 

18.9 

(30,0,10) 

8 

14568 

84% 

13.2 

18.7 

3D  1-plane 


3D  31ayer:  inplane 


3D  31ayer:  out  of 
plane 


0.012 


0.008 

0.019 


0.008 

0.022 


0.008 


Fig.  8.  The  reconstructed  absorption  coefficient  distribution  for  the  cylindrical  object  with  a 
spherical  absorption  inhomogeneity  (diameter  of  15mm  and  contrast  2:1  with  respect  to 
background)  located  at  x,  y  and  z  locations  (a)  (0,0,0),  (b)  (30,0,0)  and  (c)  (30,0,10).  The  three 
columns  of  images  show  the  results  achieved  with  the  three  different  data  collection  schemes. 
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this  discretization  and  the  associated  integration  of  the  shape  functions,  the  resolution  of  the 
mesh  and  the  associated  optical  properties  will  affect  these  results.  For  example,  if  the 
absorption  coefficient  is  much  smaller,  then  lower  resolution  meshes  may  be  adequate,  as  the 
problem  becomes  more  energy  conserving,  whereas  for  a  higher  absorption  or  scattering 
problem,  a  higher  resolution  mesh  will  be  needed  to  ensure  numerical  accuracy  within  each 
FEM  element  for  a  lossy  problem.  Note  also  that  for  spectral  reconstruction  [3]  with  six 
wavelengths  data,  each  iteration  takes  about  30  sec  for  1785  nodes  mesh. 

For  a  heterogeneous  or  homogeneous  2D  case,  number  of  useful  singular  values,  which 
are  above  the  noise  level  (1%  in  amplitude)  showed  similar  trends  and  behavior  with 
increasing  numbers  of  measurements,  as  evident  from  Fig.  3(a).  Further,  the  percentage  of 
useful  measurements  (useful  singular  values)  drops  exponentially  as  the  number  of 
measurements  is  increased,  Fig.  3(b).  It  is  worth  noting  that  for  a  heterogeneous  model,  since 
light  propagation  becomes  more  complex,  and  in  this  case  more  diffusive,  the  total  number  of 
useful  measurements  is  slightly  lower  than  that  of  homogeneous  model.  In  this  work,  useful 
singular  values  are  defined  as  the  ones  which  are  above  noise  level  (1%).  This  is  used  only  for 
optimizing  the  parameters  used  in  the  reconstruction  procedure,  but  in  the  actual 
reconstruction  procedure,  regularization  is  used  to  reduce  the  condition  number. 

Next,  the  effect  of  the  2D  mesh  resolution  was  investigated,  for  it’s  impact  upon  the 
number  of  independent  available  measurements.  From  Fig.  3(c),  it  is  evident  that  if  the 
degrees  of  freedom  (mesh  resolution)  in  the  forward  problem  is  less  than  the  total  number  of 
measurements,  then  increasing  the  number  of  measurements  does  not  increase  the  number  of 
independent  measurements  (i.e.  the  rank),  since  the  rank  is  predominantly  restricted  by  the 
number  of  nodes  in  the  mesh.  For  example,  given  a  system  from  which  only  240 
measurements  are  available,  any  mesh  which  has  a  resolution  of  240  nodes  or  more  will  give 
the  same  number  of  independent  measurements.  Therefore  no  additional  measurements  can  be 
gained  in  terms  of  independent  information  by  increasing  the  mesh  resolution.  Given  a  2D 
mesh  of  1785  nodes,  for  example,  no  considerable  gain  in  independent  data  can  be  obtained 
when  the  number  of  measurements  are  increased  above  1560  (40  source  and  detectors).  At 
this  point,  it  will  be  worth  remembering  that,  in  real  time  there  is  a  physical  constraint  on 
number  of  measurements,  because  of  the  physical  geometry  and  fiber  size.  To  take  an 
example,  for  a  circular  test  phantom  of  86  mm  diameter  and  fiber  of  6  mm  diameter,  no  more 
than  40  fibers  (which  corresponds  to  1560  measurements)  can  be  arranged  around  the  outer 
boundary  of  domain.  However  this  issue  becomes  more  important  perhaps  for  non-contact 
imaging  systems  in  which  the  number  of  source-detector  locations  can  be  arbitrarily  large. 

Using  a  2D  mesh  of  1785  nodes,  the  effect  of  an  increase  in  the  reconstruction  (regular 
pixel)  basis  resolution  upon  reconstruction  was  investigated  [Fig.  3(d)].  An  increase  in  pixel 
basis  elements  increases  the  number  of  useful  singular  values,  but  there  is  no  significant 
improvement  in  the  pixel  basis  from  30x30  (900  elements)  to  40x40  (1600  elements).  This  is 
very  interesting,  since  one  would  assume  that  fewer  degrees  of  freedom  for  the  inverse 
problem  would  produce  a  better  solution.  But  although  the  problem  may  become  better  posed, 
the  rank  will  be  similar  to  that  shown  in  Fig.  3(d).  However,  these  results  indicate  the  best 
possible  resolution  obtainable  is  by  using  the  40  x  40  pixel  basis  and  again  these  results  will 
be  dependent  on  the  physical  problem  dimension  and  level  of  complexity.  Figure  4  shows  that 
increasing  the  number  of  measurements  for  a  2D  model  increases  the  sensitivity  of  the 
problem,  as  evident  from  magnitude  plot  of  the  Jacobian  (calculated  from  1785  nodes  mesh). 
Also  shown  in  Fig.  4  is  a  normalized  plot,  relative  to  the  central  value,  and  indicates  that  for 
fewer  number  of  measurements,  the  sensitivity  is  maximal  near  the  boundary  and  lower  at  the 
center,  as  expected.  By  increasing  the  number  of  measurements,  eventually  the 
hypersensitivity  near  to  the  boundary  reduced  and  the  sensitivity  became  uniform  regardless 
of  distance  from  boundary.  Finally,  it  is  observed  that  increasing  the  number  of  measurements 
above  552  (24  sources  and  detectors)  did  not  result  in  any  further  improvement  in  the 
sensitivity  distribution. 

For  the  3D  model,  Table  2  shows  that  three  layers  of  out-of-plane  measurements  yields  a 
higher  number  of  useful  singular  values,  but  the  useful  percentage  of  the  total  measurements 
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was  below  15%.  An  increase  in  number  of  measurements  means  more  data  acquisition  time 
and  more  computation  time.  Non-linear  iterative  image  reconstruction  procedures  in  NIR 
imaging  use  repeated  calculation  of  the  forward  data.  Therefore  increasing  the  number  of 
sources  and  measurements  substantially  increases  the  computation  time.  In  comparing  the 
three  layer  in-plane  and  three  layer  out-of-plane  data  collection  strategies,  having  more  than 
three  times  the  measurements  in  the  latter  case  improves  the  number  of  useful  singular  values 
only  by  22%.  The  improvement  in  the  number  of  useful  singular  values  is  not  significant  if  the 
data  acquisition  time  is  considered  as  well  as  the  computation  time.  The  magnitude  of  the 
singular  values  indicates  the  importance  of  that  eigenvector  in  the  image  space,  which  is 
directly  related  to  reconstructed  image  contrast  that  can  be  achieved.  To  compare  the 
magnitude  of  the  largest  singular  value,  even  though  it  is  at  its  highest  for  the  three  layer  out- 
of-plane  strategy,  it  should  be  noted  that  only  3  of  the  singular  values  are  above  164 
(magnitude  of  largest  singular  value  of  3D  3 layer-in-plane),  indicating  that  there  would  not  be 
dramatic  differences  in  the  reconstructed  image  contrast  in  both  these  cases.  If  the  magnitude 
of  largest  singular  value  in  2D  and  3D  are  compared,  in  2D  the  magnitude  is  higher,  whereas 
the  number  of  useful  singular  values  are  lower  than  3D,  indicating  that  the  modes  that 
contribute  to  image  space  are  fewer  and  the  quality  of  the  reconstructed  image  in  2D  will  be 
lower  than  3D.  Even  though  magnitude  of  the  singular  values  dictate  the  contrast,  the  singular 
vectors  associated  with  it  will  tend  to  affect  the  reconstructed  image  quality.  The  magnitude  of 
the  largest  singular  value  in  the  3D  31ayer  cases  are  the  same  because  of  the  smoothness  of  the 
singular  vectors  in  the  case  of  3D  3  layer:  out-of-plane,  the  reconstructed  image  quality  is 
better  than  the  rest  cases  (Fig.  8).  The  FWHM  analysis  also  confirms  this. 

It  should  be  noted  that  there  is  always  a  trade-off  between  image  quality  and  computation 
time.  Therefore  having  out-of-plane  data  increases  the  image  resolution,  but  taking  into 
consideration  the  overall  computation  time,  this  improvement  is  perhaps  not  so  significant. 
The  computation  time  per  iteration  is  high  in  the  case  of  out-of-plane  data  (computation  time 
per  iteration:  2D  problem  -  70  sec;  single-layer  -  289  sec;  three  layer:  in-plane  -  573  sec; 
three  layer:  out-of-plane  -  1821  sec). 

Figure  5  indicates  that  for  the  3D  model  with  a  single  measurement  plane  case,  the  total 
sensitivity  is  higher  near  the  boundary,  as  compared  to  the  three  plane  data  case  and  by 
increasing  the  number  of  measurements  the  sensitivity  near  the  boundary  is  decreased.  The 
results  show  that  although  the  sensitivity  is  still  higher  at  the  boundary  with  three  planes  of 
data  acquired,  there  is  no  significant  difference  in  the  sensitivity  pattern  observed  between 
three  layer  in-plane  or  out-of-plane  strategies. 

Since  only  one  component  of  the  full  Jacobian  matrix,  J2  in  Eq.  (2),  has  been  examined 
here,  images  have  also  been  reconstructed  for  qa  using  log  amplitude  data  for  a  2D  forward 
mesh  of  1785  nodes  and  a  reconstruction  basis  30  by  30  pixel  basis.  Noisy  simulated  data 
were  generated  for  various  radial  positions  of  the  absorption  inhomogeneity  with  a  contrast  of 
2,  relative  to  the  background  and  having  a  diameter  of  10  mm.  The  log  of  RMS  error  was 
calculated  as  the  difference  in  the  original  and  the  reconstructed  horizontal  cross-sections  of 
each  image  (Fig.  6)  as  a  function  of  number  of  measurements  and  these  were  plotted  in  Fig.  7. 
The  results  show  that,  as  evident  from  Fig.  7,  although  there  is  a  decrease  in  the  RMS  error  as 
the  number  of  measurements  is  increased,  the  improvement  in  the  reconstructed  images  is  not 
significant  for  measurements  greater  than  552  (corresponding  to  24  fibers).  However,  for  a 
central  anomaly,  the  RMS  error  continued  to  decrease  with  increasing  number  of 
measurements,  whereas  for  an  anomaly  near  the  boundary  the  RMS  error  does  not  improve 
more  than  0.5%  with  respect  to  552  measurements. 

To  study  the  effect  of  data  collection  strategies  on  the  3D  reconstructed  image,  the  FWHM 
(Full  Width  at  Half  Maximum)  of  the  peaks  for  all  the  reconstructed  cases  have  been 
calculated  and  compared,  Table  3.  As  the  inhomogeneity  moves  from  the  centre  towards  the 
boundary,  the  FWHM  reduces  for  both  of  the  three  layer  cases  and  it  remains  approximately 
the  same  for  the  single  layer  case.  For  example,  when  the  inhomogeneity  is  placed  at  (30,0,0), 
Fig.  8(b),  the  FWHM  (in  the  X-cross  section)  values  for  single  layer  is  17.2mm  and  for  the 
three-layers  in-plane  and  out-of-plane  strategies  is  13.1mm  and  13.6mm  respectively.  It  is 
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evident  from  the  reconstruction  examples  that  the  quantitative  values  of  the  inhomogeneities 
increase  as  the  object  moves  from  the  centre  to  boundary,  which  is  in  close  match  with 
Jacobian  analysis  above.  Reconstruction  of  absorption  using  single  layer  data,  is  not  accurate, 
in  a  case  where  the  anomaly  is  not  presented  in  the  imaging  plane,  such  a  case  results  are 
presented  in  Fig.  8(c).  In  this  case,  single-layer  reconstructed  image  shows  the  inhomogeneity 
at  a  false  position  (reconstructed:  (30,0,0);  actual:  (30,0,10)).  Most  of  the  3D  NIR  studies 
indicate  that,  the  quantitative  accuracy  of  the  images  will  be  poor  due  the  partial  volume 
effect  in  three  dimensions[13,16,17]  and  these  quantification  can  be  greatly  improved  by  the 
use  of  more  sophisticated  regularization  and  the  addition  of  penalty  terms  into  Eq.  (3). 

5  Conclusions 

In  this  investigation,  the  mesh  resolution  and  numerical  accuracy  in  the  2D  and  3D  forward 
problems  were  examined,  using  specific  data-collection  geometries.  Several  choices  such  as 
domain  size,  optical  properties  and  anomaly  position  and  size  were  kept  fixed,  relative  to 
typical  breast  cancer  imaging  situations.  It  was  shown  that  increasing  the  number  of 
measurements  increases  the  total  amount  of  information  available,  and  these  specifically 
enhance  the  recovery  of  the  central  region  of  the  model,  regardless  of  dimensionality.  Further, 
by  increasing  the  number  of  measurements,  the  rank  of  the  problem  (i.e.  amount  of 
independent  useful  information)  may  not  increase  if  the  degrees  of  freedom  (i.e.  number  of 
nodes  in  the  mesh)  are  low.  Reconstruction  basis  plays  an  important  role  in  the  inverse 
problem  and  it  has  been  found  that  a  pixel  basis  of  30  x  30  is  optimal  for  a  typical  breast 
imaging  problem. 

More  specifically  for  a  3D  imaging  problem,  this  work  has  shown  the  benefits  and 
drawbacks  of  multi-plane  data  collection  as  well  as  the  use  of  in-plane  versus  out-of-plane 
data  measurements  strategies.  It  has  been  shown  that  the  use  of  single-plane  of  data  in  a  3D 
model  is  perhaps  adequate,  in  terms  of  image  quality,  computation  time  and  data  collection 
time,  if  the  anomaly  being  imaged  is  within  the  plane  of  measurements.  However,  if  prior 
information  such  as  plane  of  interest  is  not  known,  it  has  been  shown  that  multi-plane  data  is 
crucial.  The  use  of  in-plane  and  out-of-plane  data  has  been  addressed  and  is  shown  that 
although  the  use  of  out-of-plane  data  provides  more  independent  and  useful  information  for 
image  reconstruction,  the  magnitude  of  this  additional  information  does  not  provide  enough 
advantages  worth  the  data  acquisition  and  image  computation  time. 

Finally  it  is  worth  noting  that  the  3D  study  has  been  limited  to  16  source/detection  fibers 
per  plane.  The  addition  of  more  measurement  fibers  and/or  investigation  of  a  different  image 
reconstruction  basis,  such  as  those  performed  for  the  2D  problem  can  be  easily  extended  for 
the  presented  3D  problem.  The  technique  and  analysis  described  here  can  be  used  as  a  tool  to 
improve  resolution  and  contrast,  given  prior  information  about  the  domain  being  imaged.  This 
specific  study  was  undertaken  to  better  understand  the  parameters  and  capabilities  of  existing 
breast  imaging  system  at  Dartmouth  and  to  focus  on  software  improvements  which  may 
increase  its  recovery  of  lesion  information. 
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Diffuse  optical  tomography  (DOT)  involves  estimation  of  tissue  optical  properties  using  noninva- 
sive  boundary  measurements.  The  image  reconstruction  procedure  is  a  nonlinear,  ill-posed,  and 
ill-determined  problem,  so  overcoming  these  difficulties  requires  regularization  of  the  solution. 

While  the  methods  developed  for  solving  the  DOT  image  reconstruction  procedure  have  a  long 
history,  there  is  less  direct  evidence  on  the  optimal  regularization  methods,  or  exploring  a  common 
theoretical  framework  for  techniques  which  uses  least-squares  (LS)  minimization.  A  generalized 
least-squares  (GLS)  method  is  discussed  here,  which  takes  into  account  the  variances  and  covari¬ 
ances  among  the  individual  data  points  and  optical  properties  in  the  image  into  a  structured  weight 
matrix.  It  is  shown  that  most  of  the  least-squares  techniques  applied  in  DOT  can  be  considered  as 
special  cases  of  this  more  generalized  LS  approach.  The  performance  of  three  minimization  tech¬ 
niques  using  the  same  implementation  scheme  is  compared  using  test  problems  with  increasing 
noise  level  and  increasing  complexity  within  the  imaging  field.  Techniques  that  use  spatial-prior 
information  as  constraints  can  be  also  incorporated  into  the  GLS  formalism.  It  is  also  illustrated  that 
inclusion  of  spatial  priors  reduces  the  image  error  by  at  least  a  factor  of  2.  The  improvement  of 
GLS  minimization  is  even  more  apparent  when  the  noise  level  in  the  data  is  high  (as  high  as  10%), 
indicating  that  the  benefits  of  this  approach  are  important  for  reconstruction  of  data  in  a  routine 
setting  where  the  data  variance  can  be  known  based  upon  the  signal  to  noise  properties  of  the 
instruments.  ©  2007  American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.2733803] 

Key  words:  near  infrared,  diffuse  optical  tomography,  inverse  problems,  least-squares  minimiza¬ 
tion 


I.  INTRODUCTION 

Image  reconstruction  methods  used  in  diffuse  optical  tomog¬ 
raphy  (DOT)  are  mainly  dependent  on  the  type  of  data,  the 
diffuse  light  model,  and  the  number  of  available  anatomical/ 
spectral  priors.  There  are  numerous  reconstruction  tech¬ 
niques  available  in  the  literature  depending  on  the 
application. 1-5  Yet  despite  the  volume  of  work  in  this  area 
there  is  no  single  investigation  with  a  direct  comparison  of 
the  least-squares  (LS)  minimization  techniques  using  the 
same  implementation  scheme,  especially  in  terms  of  data 
noise  level  and  complexity  in  the  test  fields.  Most  of  the 
comparisons  in  the  literature  have  been  in  terms  of  imple¬ 
mentation  of  minimization  and  convergence  rates  of  one  or 
two  techniques  at  hand. 1-5  This  work  addresses  this  problem 
and  compares  minimization  methods  (more  specifically  dif¬ 
ferent  types  of  regularization)  with  the  same  implementation 
scheme  for  a  direct  quantitative  comparison.  Moreover,  us¬ 
age  of  weight  matrices  in  the  regularization  which  include 
the  variance  and  covariance  properties  of  data  and  image 
space  are  extensively  explored  here.  A  new  covariance  form 


borrowed  from  meteorological  studies  is  introduced  and 
proven  to  be  effective  for  reconstructing  highly  noisy  data  in 
the  generalized  theoretical  frame  work. 

Near  infrared  DOT  involves  reconstructing  images  of  op¬ 
tical  properties  from  transmission  measurements  using  wave¬ 
lengths  from  650-1000  nm  to  interrogate  tissue. 1,6-8  Optical 
absorption  and  scattering  images  obtained  using  multiple 
wavelengths  can  be  used  to  estimate  tissue  hemoglobin,  wa¬ 
ter  concentration,  scattering  amplitude,  and  scattering 
power.8  To  overcome  the  inherent  low-spatial  resolution  in 
DOT,  there  is  a  considerable  interest  in  developing  hybrid 
systems,9-27  which  use  the  spatial  mapping  of  one  system  as 
the  template  for  DOT.  Image  formation  from  the  data  col¬ 
lected  by  these  (stand-alone/hybrid)  systems  involves  solv¬ 
ing  an  inversion  problem.  This  article  describes  LS  minimi¬ 
zation  techniques  to  solve  the  inverse  problem  and  to 
quantitatively  compare  their  performance  in  a  systematic  se¬ 
ries  of  simulations.  The  inverse  problem  (image  reconstruc¬ 
tion  procedure)  in  DOT  is  known  to  be  a  nonlinear,  ill-posed, 
and  ill-determined  problem,2  and  to  solve  such  a  problem,  a 
regularization  term  must  be  added  to  constrain  the  solution 
space  in  order  to  obtain  a  meaningful  image.  There  are  many 
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Forward  Problem  Inverse  Problem 


Fig.  1 .  An  illustration  of  the  forward  and  inverse  problem  in  diffuse  optical 
tomography  is  shown  (see  Ref.  64),  where  (a)  the  data  y  is  estimated  given 
values  of  /jia  and  /jl's  and  source/detector  positions.  In  the  inverse  problem 
(b),  the  values  of  i±a  and  jm's  must  be  obtained  given  a  set  of  measurements 

(y). 


regularization  methods  available  in  the  literature  and  this 
work  focuses  on  the  fact  that  most  LS  techniques  presented 
in  the  literature  can  be  encompassed  within  a  generalized 
theoretical  framework,  which  includes  a  regularization  ma¬ 
trix  that  is  based  upon  weights  from  the  data  and  parameter 
variances.  Note  that  Appendix  A  gives  the  terminology  used 
in  this  work  along  with  definitions  of  symbols. 

Because  of  the  interest  in  using  spatial  information  de¬ 
rived  from  conventional  imaging  modalities  in  the  DOT  in¬ 
verse  problem,  a  number  of  methods  have  been  presented  in 
the  literature.9-27  These  techniques  were  initially  proposed  by 
Barbour  et  al.  and  Schweiger  et  al.  and  used  in  to  improve 
the  quantitative  outcome  of  reconstructed  images.  Ntziach- 
ristos  et  al. 14  used  the  magnetic  resonance  (MR)  information 
to  divide  the  imaging  domain  into  tumor  and  nontumor  re¬ 
gions  to  make  the  problem  better  posed.  Li  et  al.11  used  an 
x-ray  tomosynthesis  volume  to  segment  the  breast  into  dif¬ 
ferent  subregions  and  used  different  regularization  param¬ 
eters  depending  on  the  size  of  the  subregions.  Recently  Gu- 
ven  et  al.24  proposed  a  Bayesian  frame  work  to  include 
spatial  prior  information  in  an  effective  way  which  will  not 
bias  the  image  reconstruction  problem  to  imperfect  anatomi- 
cal  priors.  Pogue  and  Paulsen,  Brooksby  et  al.  ’  ’  and 
Yalavarthy  et  al.26  have  extended  these  approaches  for  the 
use  of  anatomical  prior  information  in  which,  depending  on 
the  connectivity  and  size  of  the  subregion,  the  regularization 
term  was  scaled.  Even  though  the  effect  of  imperfect  spatial 
prior  information  on  the  image  reconstruction  is  a  very  active 
research  area,  ’  ’  it  was  assumed  here  that  the  spatial  pri¬ 
ors  were  perfect.  Other  ongoing  studies  are  examining  this 
more  complex  issue. 


II.  DOT  FORWARD  PROBLEM 

DOT  involves  solving  a  model  (forward)  and  estimation 
(inverse)  problem,  sequentially  as  illustrated  in  Fig.  1.  In  this 
section,  the  forward  problem  is  described,  which  involves 
generating  the  measurement  data,  for  a  given  set  of  optical 
property  estimates  within  the  tissue,  using  a  finite  element 
solution  to  the  diffuse  transport  equation. 


Light  propagation  in  a  turbid  elastic- scattering  media,  like 
tissue,  is  treated  as  “neutral-particle  transport”  rather  than 

“wave  propagation”  and  in  the  frequency  domain,  the  diffu- 

2  28 

sion  equation  is  used,  which  is  given  by  ’ 

-  V  .  D(r)  V  T>(r,  o>)  +  [jma( r)  +  zWc]<l>( r,o>)  =  Q0( r,  w), 

(i) 


where  <L(r,co)  is  the  photon  density  at  position  r  and  the 
light  modulation  frequency  is  given  by  co  (<u=27t/,  in  this 
work  /=  100  MHz).  The  isotropic  source  term  is  represented 
by  Q0( r ,  a) )  and  the  speed  of  light  in  tissue  by  c,  which  is 
constant  here.  jma( r)  is  the  optical  absorption  coefficient  and 
D(r)  is  the  optical  diffusion  coefficient,  which  is  defined  as 


D(  r)  = 


1 

3[juLa(r)  +  jmfs(r)y 


(2) 


where  /z^(r)  is  the  reduced  scattering  coefficient,  which  is 
defined  as  jul's  =  /z^(l  -g).  jms  is  the  scattering  coefficient  and  g 
is  the  anisotropy  factor.  A  Robin  (type-III)  boundary  condi¬ 
tion  is  applied  to  model  the  refractive-index  mismatch  at  the 
boundary.  The  measured  data  for  a  frequency  domain  sys¬ 
tem  are  the  amplitude  and  phase  of  the  transmitted  signal.  If 
F  is  the  forward  model  [finite  element  method  (FEM)  in 
here]  which  gives  the  fluence  at  every  point,  then  the 
modeled  data  G(/z)  can  be  obtained  by  sampling  the 
forward  model  at  the  boundary  given  internal  spatial  distri¬ 
butions  of  optical  properties  and  source-detector  locations, 
where  /z  represents  the  parameters  {/z=[D(r)  ;/za(r)]}, 

G(/ul)  =  S{F(jul)}.  (3) 


The  details  of  the  FEM  formulation  of  the  forward  model 
are  given  in  Refs.  30-32.  The  results  presented  are  restricted 
to  frequency-domain  data,  more  specifically  data  (y)  is  the 
natural  logarithm  of  the  amplitude  (A)  and  phase  (6)  of 
the  frequency-domain  signal.  Defining  A  and  6  in  terms 
of  modeled  data,  A  =  ^Re{G(/z)}2+Im{G(/z)}2  and  0 
=  tan_1(Im{G(/z)}/Re{G(/z)}).  The  Jacobian  (/),  which  gives 
the  rate  of  change  of  modeled  data  with  respect  to  param- 
eters,  is  calculated  using  the  adjoint  method.  Even  though 
the  actual  parameters  being  estimated  are  D( r)  and  jma( r),  the 
results  are  presented  in  terms  of  jna( r)  and  jul's( r),  which  are 
spectroscopically  more  meaningful. 


III.  LEAST-SQUARES  MINIMIZATION  TECHNIQUES 

This  section  outlines  several  different  minimization 
schemes  used  in  this  work.  These  techniques  are  used  to 
solve  the  inverse  problem  [Fig.  1(b)],  which  is  achieved  by 
minimizing  the  objective  function  (12)  over  the  range  of  /z. 
Minimizing  the  objective  function  can  be  achieved  by  sev¬ 
eral  different  approaches.  The  most  common  approaches  in¬ 
volve  obtaining  repeated  solutions  of  the  forward  model  and 
recomputation  of  the  Jacobian  (/)  (and  its  inversion)  at  every 
iteration  because  of  the  nonlinear  nature  of  the  problem. 
There  are  also  gradient-based  optimization  schemes  avail- 
able  in  the  literature  ’  to  minimize  the  objective  function 
which  does  not  require  an  explicit  inversion  of  the  Hessian 
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matrix.  In  this  work  direct  methods,  known  as  full-Newton 
approaches,  are  employed  in  minimization  for  all  the  regu¬ 
larization  techniques  used  for  a  fair  comparison.  LS  minimi¬ 
zation  has  the  effect  of  reducing  high  frequency  noise,  lead¬ 
ing  to  smooth  images  of  optical  properties.  Total  variation 
methods  and  variants  of  this  are  used  to  obtain  edge  preser¬ 
vation  in  reconstructed  images.27,35  Solving  the  inverse  prob¬ 
lem  using  LS  minimization  can  also  be  seen  from  a  Bayesian 
prospective  to  obtain  maximum  a  posteriori  estimate.24,36,37 
A  correlation  between  the  Bayesian  frame  work  and  LS 
minimization  techniques  is  given  in  Refs.  12,  38,  and  39,  but 
usage  of  the  Bayesian  frame  work  requires  one  to  choose  a 
particular  noise  model  for  both  data  and  image  space,  which 
might  not  reflect  the  actual  noise  characteristics  unless  some 
prior  information  is  available.  Here,  the  emphasis  is  on  LS 
minimization  techniques  with  a  focus  on  what  the  value  of 
the  regularization  method  can  be.  The  LS  methods  are  di¬ 
vided  into  two  groups:  (1)  Without  spatial  priors  and  (2)  with 
spatial  priors. 


A.  Without  spatial  priors 
1 .  Levenberg-Marquardt  minimization 

5  39 

This  approach  is  also  known  as  a  trust-region  method  ’ 
where  experimental  data  is  matched  with  modeled  data 
iteratively.40,41  The  objective  function  for  the  DOT  problem 
is  defined  as 

fi  =  {|b'-G(/u,)||2},  (4) 

where  y  is  the  data  and  G(/x)  is  the  modeled  data.  This  equa¬ 
tion  is  minimized  by  setting  the  first-order  derivative  equal  to 
zero. 

a.  First-order  condition.  Minimizing  12  with  respect  to 
/x,  which  is  achieved  by  setting  <?I2/<?/x=0, 

dfl  T 

—  =JTS=  0,  (5) 

dp, 

where  8  is  the  data-model  misfit,  £=y-G(/x),  J  is  the  Jaco¬ 
bian,  and  T  represents  the  matrix  transpose  operator. 

b.  Iterative  update  equation.  Imagine  a  sequence  of  ap¬ 
proximations  to  /x  represented  by  then  using  Taylor  series 


on  G(/X;)  and  expanding  around  /x^  gives 

G(ni)  =  G(,u.(_i)  +  JA^  +  ...  ,  (6) 

where  A /xz- = /xz  —  ytxz_  x .  Rewriting  8  utilizing  the  first  two 
terms  of  Eq.  (6)  (ignoring  the  rest,  equivalently  linearizing 
the  problem)  gives 

Sj  =  y  -  G(fii)  =y-  G(/i,_1)  -  JA/Xj  =  S M  -JA/xr  (7) 
Rewriting  Eq.  (5)  for  the  ith  iteration 
JT8t  =  0.  (8) 

Substituting  Eq.  (7)  into  Eq.  (8)  gives 

~JAfJ.t)  =  0.  (9) 

Further  simplification  leads  to  the  update  equation 
[JTJ]AjuLi  =  JT8i_l.  (10) 


When  JTJ  is  ill-conditioned,  a  diagonal  term  is  added  to 
stabilize  the  problem.  In  this  case,  the  update  equation  be¬ 
comes: 

[  JrJ  +  aljAfii  =  Jt8i_i  ,  (11) 

where  A/x,  is  the  update  for  the  parameter  in  the  ith  step. 
Note  that  a  monotonically  decreases  with  iterations  (always 
>0),  and  also  that  o^||<5||2.  The  iterative  method  (or  its 
modified  version)  is  the  commonly  used  minimization  tech¬ 
nique  in  DOT.  It  can  be  seen  from  Eqs.  (10)  and  (11),  when 
a  becomes  zero  in  Eq.  (11)  it  becomes  Eq.  (10).  It  is  also 
important  to  note  that  JTJ  is  always  symmetric,  because 
{JTJ)T=JT(JT)T=JTJ.  The  advantage  of  using  this  method  is 
in  the  simple  choice  of  a  regularization  parameter  ( a ).  The 
limitations41  of  this  method  include: 

•  JTJ  must  be  positive  definite. 

•  The  initial  guess  (/x0)  should  be  close  to  the  actual  so¬ 
lution. 

•  The  update  equation  [Eq.  (11)]  does  not  solve  the  first- 
order  conditions  unless  a= 0. 

•  Since  parameters  are  not  involved  in  the  minimization 
scheme,  the  inverse  problem  may  be  unstable. 

Even  though  JTJ  is  not  positive  definite  in  DOT,  the 
Levenberg-Marquardt  (LM)  approach  (or  its  modified  ver¬ 
sion)  has  been  used  successfully  in  a  number  of 

instances.2,6,7,28,42 


2 .  Tikhonov  minimization 

The  generalized  objective  function43,44  in  the  Tikhonov 
case  includes  parameters  in  the  minimization  function,  which 
is  defined  as 

a  =  {||y-G(M)||2  +  X||L(M-^o)l|2},  (12) 


where  X  is  the  Tikhonov  regularization  parameter  and  L  is  a 
dimensionless  regularization  matrix  (in  this  work).  Here,  /x0 
is  the  prior  estimate  of  the  optical  properties,  which  in  DOT 
has  typically  been  obtained  from  calibrating  the  data.45,46 

a.  Choice  ofX.  Rewriting  Eq.  (12),  normalizing  both 
terms  by  their  variances  yields 


^  _  I  l|y  -  G(/x)||2  +  II L(fi  -  At0)||2 


K)2 


(<W0) 


(13) 


where  cry  is  the  standard  deviation  in  the  data  y  and  a is 
the  standard  deviation  in  the  optical  properties  (or  deviation 
from  the  prior  estimate  of  optical  properties).  Note  that  the 
variance  of  data-model  misfit  [ 8=y-G(p )]  is  assumed  from 
the  data,  i.e.,  (o's)2=(cry)2-\-(crG^/ULy)2  with  (crG^/ULy)2=0  because 
synthetic  data  was  used.  Multiplying  Eq.  (13)  by  and 
comparing  the  result  with  Eq.  (12)  leads  to 


X  = 


K)2 


(14) 


which  shows  that  the  Tikhonov  regularization  parameter  (X) 
should  be  equal  to  the  square  of  the  ratio  of  the  standard 
deviation  in  data  to  the  standard  deviation  of  the  parameters. 
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This  is  a  subtle  yet  important  point,  especially  since  this 
parameter  is  rarely  defined  this  way,  and  is  most  commonly 
derived  empirically. 

b.  First-order  condition.  Minimizing  H  with  respect  to 
/z,  which  is  achieved  by  setting  dFU  dp- 0, 

dFl 

—  =  JT8-\LTL{fjL-^)  =  0.  (15) 

d/JL 

c.  Update  equation.  Rewriting  Eq.  (15)  for  the  ith  itera¬ 
tion  leads  to 

JTSi-\LTL(vi-v0)  =  0.  (16) 

Substituting  Eq.  (7)  into  Eq.  (16)  results  in 

JT(Si-i  -  J^i)  -  M/L(m,-i  +  Am,  -  Mo)  =  o.  (17) 

Further  simplification  leads  to  the  iterative  update  equa¬ 
tion 

[JrJ  +  \I/L]Am,  =  Jr<5,-i  -  M/L(m,-i  -  Mo)  •  (18) 

Note  that  LrL  is  symmetric.  The  constraint  on  the  choice 
of  L  is  that  it  must  be  positive  definite.44  In  the  absence  of 
spatial  priors,  a  common  choice  for  the  form  of  L  is  the 
identity  matrix  (/),  which  leads  to  the  update  equation 

[JrJ  +  \I]Am,  =  Jr<S,-i  -  Mm,-i  -  Vo)  ■  (19) 

Refer  to  Appendix  B  for  an  analysis  of  the  Tikhonov  regu¬ 
larization  in  terms  of  singular  values.  This  regularization 
method  is  particularly  common  for  ill-posed  problems.  The 
advantage  of  the  method,  is  that  it  includes  parameters 
within  the  minimization  scheme  which  can  be  selected  to 
improve  the  stability  of  the  solution.  Its  limitations  are  that: 

•  it  requires  a  prior  opinion  about  the  noise  characteristics 
of  the  parameter  and  data  spaces  (for  X)  and 

•  it  does  not  take  into  account  the  individual  variances  of 
the  data  points/parameters,  nor  their  covariances. 

However,  the  simplicity  of  the  approach  makes  it  attrac¬ 
tive  for  use  in  ill-posed  problems.  When  the  dynamic  range 
of  the  data  is  large  (as  in  DOT),  incorporation  of  the  maxi¬ 
mum  variance  in  the  data  will  cause  the  minimization  to  bias 
the  solution  to  specific  data  points  (e.g.  near  the  boundaries 
at  source-detector  locations  in  DOT).  To  reduce  the  effect  of 
bias,  one  can  employ  a  generalized  least  squares  (GLS)  mini¬ 
mization  scheme,  described  in  the  next  section. 


3 .  GLS  minimization 

Generalized  least  squares  minimization  schemes  have 
been  proposed  in  the  context  of  Tikhonov  minimization  in 
the  literature,  ’  ’  in  which  there  is  some  ambiguity  in  choos¬ 
ing  the  regularization  parameter  (X).  In  here,  a  direct  inclu¬ 
sion  of  weight  matrices  (which  are  inverses  of  covariance 
matrices)  in  the  minimization  scheme  was  employed  to  ex¬ 
plicitly  remove  the  dependence  of  reconstructed  image  qual¬ 
ity  on  the  choice  of  regularization  parameter.  This  type  of 
choice  leads  to  an  objective  function47,48 


a  =  {[y-G(M)FWiy-G(M)] 

+  (20) 

where  Ws  is  the  weight  matrix  for  data-model  misfit  (S)  with 
W^=[cov(^)]_1  (Appendix  A-4  of  Ref.  47).  is  the 

weight  matrix  for  optical  properties  (p-p0)  with 
=  [cov(/z-/z0)]_1  (Appendix  A-4  of  Ref.  47).  Explicit  forms 
for  these  weight  matrices  are  discussed  later.  Since  both  are 
inverses  of  covariance  matrices,  they  are  symmetric  and 
positive  definite. 

a.  First-order  condition.  Minimizing  D  with  respect  to 
/z,  which  is  achieved  by  setting  dFll  dp- 0  produces 

dFl 

—  =JTWs8-Wfl_tlo(^-fi0)  =  0.  (21) 

b.  Update  equation.  Similar  to  the  Tikhonov  approach, 
linearizing  the  problem  leads  to  the  iterative  update 

..  48 

equation 

[J,W(J  +  W^0]Am,  =  J'WmVi  -  W^Mi-1  -  Mo)  • 

(22) 


4.  Choice  of  Ws 

Since  simulated  data  were  used  here,  in  the  formation  of 
the  weight  matrix  (covariance  matrix),  it  was  assumed  that 
the  cov(  8)  is  due  to  measurement  error  only,  which  yields47 

Ws=  [cov(<?)]_1  =  {cov[y  -  G(m)]}"1  =  [cov(y)]-1,  (23) 


where  cov  represents  the  covariance  operator.  In  the  simula¬ 
tion,  typically  one  generates  the  forward  data  and  adds  noise 
to  it  to  form  synthetic  data 

y  =  G(fi)  +  oyi 7,  (24) 


where  77  is  a  random  number  vector.  Typically,  a  random 
number  generator  which  follows  a  normal  distribution  with 
zero  mean  and  unity  variance  is  used.  Here,  cry  is  the  stan¬ 
dard  deviation  of  the  data,  assuming  the  noise  is  totally  un¬ 
correlated  (white  noise)  in  which  case,  the  covariance  matrix 
becomes47 


[cov(y)]0-  = 


1 0  if  i  A  j 
1  («■>)?  if  '=./• 


(25) 


Since  synthetic  data  were  used  in  this  article,  the  weight 
matrix  for  the  data  (Ws)  becomes  diagonal.  In  the  experi¬ 
mental  case,  one  needs  to  collect  an  ensemble  of  data  sets 
from  which  a  covariance  matrix  can  be  computed.  In  this 
case,  “A”  data  sets  need  to  be  collected  using  the  same  phan¬ 
tom  (different  homogeneous  phantoms  need  to  be  used  for 
different  signal  levels),  where  N  needs  to  be  a  large  number. 
From  this  ensemble  of  {y}, 


{y}=y  +  (y}> 


(26) 


where  y  is  the  true  or  mean  value  of  data  and  {v}  is  pertur- 
bation  due  to  noise.  This  leads  to 
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[cov(y)]  =  [cov(j)]  =  {y}{y}r=  —  ^  (27) 

substituting  Eq.  (27)  in  Eq.  (23)  gives  Ws.  Note  that  in  Eq. 
(23),  it  was  assumed  the  cov(<5)  is  due  to  measurement  error, 
which  is  also  true  in  the  case  of  experimental  data,  as  the 
data  are  calibrated  to  remove  the  offset  and  match  the  mod- 
eled  data.45’46 


5.  Choices  of 

Here,  two  forms  were  considered  to  highlight  the  versa¬ 
tility  of  the  procedure,  even  though  many  other  forms  of  the 
covariance  matrix  can  exist. 

a.  Analytical  covariance  form.  Borrowed  from  the  me¬ 
teorological  studies,  assuming  the  parameter  field  obeys  the 
Helmholtz  equation,  an  analytical  form  (for  one-dimensional 
infinite  domain  case)  for  the  covariance  matrix  is47 

[cov(/t  -  Ho)lj  =  «W2( 1  +  r~i]e~rif 7’  (28) 


where  r ^  is  the  separation  distance  between  locations  and  /  is 
the  correlation  length  scale.  (cr^^J2  is  the  expected  variance 
for  jul-/i0.  In  this  case,  the  weight  matrix  is  constructed  from 
^-Mo=[cov(/i-^0)]-1 

b.  Local  Laplacian  form.  Here,  the  weight  matrix  is 
formed  directly  using  a  local  Laplacian  operator5,49,50  be¬ 
tween  neighboring  locations,  where 

W^(=[\/(a^o)2]MTM, 


where  M  is  the  local  Laplacian  matrix,  which  is  defined  as 


o 

- 1 


2  Mu 


if  i  and  j  are  not  neighbors 
if  i  and  j  are  neighbors 

if  i=j 


(29) 


where  i  and  j  represent  the  node  numbers  of  the  FEM  mesh, 
which  in  turn  become  the  indices  of  the  local  Laplacian  ma¬ 
trix  (M).  The  diagonal  terms  in  M  gives  the  total  number  of 
immediately  connected  nodes. 

Computation  of  requires  an  estimate  of  variance  of 

parameters  [(er^^)2],  as  is  the  case  for  calculation  of  the 
Tikhonov  regularization  parameter  [Eq.  (14)].  The  expected 
variance  can  be  computed  in  many  ways:  the  most  common 
method  for  imaging  problems  is  estimation  from  the  litera¬ 
ture.  For  example,  the  optical  contrast  between  tumor  and 
normal  breast  tissue  is  around  50%-400%  (Refs.  51  and  52), 
which  gives  the  expected  standard  deviation  (o^^)  in  the 
optical  properties,  and  can  be  used  to  compute  variance.  The 
calibration  of  the  experimental  data  is  capable  of  giving  a 
very  good  estimate  of  normal  tissue  optical  properties.45,46 
Note  that  weight  matrix  containing  the  expected  variance 
will  not  impose  a  hard  constrain  on  the  expected  optical 
properties,  but  discourages  update  values  (A/z)  which  are 
above  these  expected  deviation  in  a  given  iteration. 

The  advantages  of  the  GLS  approach  are  that: 


•  It  accounts  for  covariance  among  the  parameters  and 
data  points. 

•  It  also  allows  the  individual  data  points/parameters  to 
have  different  noise  characteristics  (variances). 

•  It  constrains  the  problem  through  the  weight  matrices, 
to  produce  stable  solutions. 

The  limitations  of  the  procedure  are: 

•  It  requires  prior  knowledge  about  the  noise  characteris¬ 
tics  of  parameter  and  data  spaces. 

•  The  weight  matrices  may  necessitate  computation  of  the 
inverse  of  covariance  matrices  (increasing  run  time  and 
memory  requirements). 

•  It  can  generate  unstable  solutions  when  unreasonable 
constraints  are  inadvertently  applied. 


B.  With  spatial  priors 

Overall,  the  LS  minimization  schemes  using  spatial  priors 
can  be  broadly  classified  into  two  approaches:  (1)  soft-priors 
and  (2)  hard-priors.  The  following  two  subsections  will  dis¬ 
cuss  these  two  approaches. 


1 .  Soft-priors 


In  this  approach,  the  regularization  matrix  L  in  the 
Tikhonov  approach  [Eq.  (18)]  encodes  the  spatial 
information.21,26  Previous  results  have  shown  that  using  the 
spatial  priors  in  this  fashion  do  not  bias  the  image  estimate 
when  the  prior  information  is  imperfect.26  Typically,  the  con¬ 
ventional  image  is  segmented  into  different  regions  depend¬ 
ing  on  tissue  type  to  generate  the  spatial  constraints.  The  L 
matrix  relates  each  nodal  optical  property  in  the  numerical 
model  to  the  other  nodes  in  that  region.26  Two  possible  forms 
are  indicated  later. 
a.  Laplacian  form11 


h= 


0 

I  -  UN 
1 


if  i  and  j  are  not  in  the  same  region 
if  i  and  j  are  in  the  same  region 
if  i=j 


(30) 


where  N  is  the  number  of  sampling  points  (e.g.,  nodes  in  a 
FEM  mesh)  in  that  region. 
b.  Helmholtz  form16 


V 


o 

-  1  l[N+(Kh)2] 
1 


if  i  and  j  are  not  in  the  same  region 
if  i  and  j  are  in  the  same  region 
if  i=j 


(31) 


where  N  is  the  number  of  nodes  in  that  region,  k-\U  with  / 
being  the  covariance  length  and  h  is  the  distance  between  the 
nodes. 


2 .  Hard-priors 

In  the  hard-prior  approach,  also  known  as  a  parameter- 
reduction  technique,  the  number  of  parameters  to  be  esti- 
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Fig.  2.  The  chosen  optical  property  distribution/domain  for  the  generation 
of  synthetic  data  is  shown.  The  diameter  of  the  domain  was  86  mm. 


mated  becomes  the  number  of  regions  segmented  from  the 
other  imaging  modality  (spatial  priors).26,53  Even  though  the 
number  of  parameters  to  be  estimated  reduces  considerably 
(relative  to  soft-priors),  the  problem  can  still  be  ill-posed,2  so 
a  LM  approach  was  used  [Eq.  (11)]  in  this  case  due  to  its 
simplicity.  The  main  advantages  of  the  method  are: 

•  The  problem  is  overdetermined,  which  also  implies  JTJ 
is  positive  definite. 

•  It  is  computationally  efficient. 

The  limitations  include: 

•  The  effect  of  error  or  uncertainty  in  the  spatial  priors 
can  be  amplified  by  the  technique. 

•  The  DOT  problem  may  still  be  ill-posed  (and  ill- 
conditioned)  after  the  constraints  are  added.2 


3 .  Important  notes  about  minimization  schemes 

There  are  additional  important  points  about  these  minimi¬ 
zation  schemes. 

•  The  weight  matrices  (Ws  and  )  in  the  GLS 

scheme  are  computed  before  the  iterative  reconstruction 
procedure  begins  and  are  invariant  during  the  iterative 
process.  The  same  is  true  of  the  soft-priors  L-matrix 
calculations. 

•  The  first-order  conditions  [Eqs.  (5),  (15),  and  (21)]  de¬ 
rived  by  minimizing  the  objective  functions  [Eqs.  (4), 
(12),  and  (20)]  in  all  minimization  schemes  appear  on 
the  right-hand  side  (rhs)  of  the  update  equations  [Eqs. 
(11),  (18),  and  (22)]  which  means  that  only  when  the 
rhs,  has  reached  zero,  the  solution  reached  the  global 
minimum. 

•  Computation  of  weight  matrices,  L  matrices  and  the 
Tikhonov  regularization  parameter,  requires  a  prior 
opinion  about  the  variances  of  the  parameters  and  data. 
Here,  only  the  best  prior  estimates  are  used,  which 
means  that  the  actual  variances  of  the  parameter  and 
data  spaces  are  used  in  the  reconstruction  procedure. 
Variation  from  the  best  prior  values  can  be  examined 
also,  to  observe  the  effect  of  priors,  but  that  work  is 
beyond  the  scope  of  the  present  article. 

•  When  spatial  priors  are  used  in  this  study  (as  well  as  in 
most  studies),  it  is  assumed  that  they  are  perfect.  The 


Actual  LM  Tikhonov  GLS-AC  GLS-LL 

f 


0.008  0.01  0.012  0.014  0.016  0.018  0.02  0.022 

•  •  •  •  • 

(a)  1  1.5  2  2.5  3 

Laplacian  Helmholtz  Hard-Priors 

CC9 


0.008  0.01  0.012  0.014  0.016  0.018  0.02  0.022 


(b)  1  1.5  2  2.5  3 


Fig.  3.  Reconstruction  results  (top  of  the  first  row,  abbreviations  are  given 
in  Appendix  A)  are  shown  using  noiseless  data  (bias  calculations)  (a)  with¬ 
out  spatial  priors  and  (b)  with  spatial  priors.  The  top  row  contains  images  of 
yu,fl  and  bottom  row  shows  jul's  images. 

effect  of  spatial  prior  uncertainty  on  the  DOT  inverse 
problem  is  discussed  in  Refs.  23,  24,  and  26  and  is  the 
subject  of  ongoing  study. 

•  The  covariance  lengths  associated  in  the  weight  matrix 
[GLS-analytical  covariance  (AC)  form,  Eq.  (28)],  and 
the  L  matrix  [Helmholtz  form,  Eq.  (31)]  calculations  are 
chosen  to  be  10  and  5  mm,  respectively.  The  effect  of 
covariance  length  on  the  image  reconstruction  is  dis¬ 
cussed  in  Ref.  26. 

•  In  the  LM  approach  [Eq.  (11)],  the  Jacobian  is  normal¬ 
ized  by  its  optical  properties.  Also  a  was  chosen  ini¬ 
tially  to  be  1  and  it  was  reduced  by  a  factor  of  10025  at 
every  iteration  and  multiplied  by  the  maximum  of  the 
diagonal  values  of  JTJ.  The  normalization  procedure  is 
described  in  Ref.  54.  Moreover,  eight  iterations  were 
chosen  for  all  the  LM  reconstructions,  as  it  has  been 
shown  in  the  literature  that  after  this  iteration,  error  in 
the  optical  properties  increases  for  this  particular  prob¬ 
lem  and  algorithm.55,56  This  inherent  instability  can  be 
attributed  to  the  fact  that  JTJ  is  not  positive  definite  in 
DOT. 

•  For  simplicity,  all  the  reconstruction  algorithms  are 
tested  only  in  the  two-dimensional  case.  Comparison  of 
three-dimensional  reconstructions  are  left  for  future 
investigations. 

4.  Special  cases  of  GLS  minimization 

The  update  equation  for  the  GLS  scheme,  Eq.  (22),  turns 
into  the  Tikhonov  case  [Eq.  (18)]  when  WS=I  and  W 
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Table  I.  Mean  and  standard  deviation  of  the  reconstructed:  (a)  /jia  and  (b)  jul's  values  (in  mm-1)  for  different 
regions  [labeled  in  first  column  of  Fig.  3(a)]  recovered  with  data  having  0%,  5%,  10%  noise  for  images  shown 
in  Figs.  3-5. 


Noise 

Method 

level 

Region-0 

Region- 1 

Region-2 

Actual 

— 

0.01 

0.02 

0.01 

LM 

0% 

0.0101±0.001 

0.0172±0.0023 

0.0105±0.0005 

5% 

0.0102±0.0016 

0.0125±0.0016 

0.0123±0.0011 

10% 

0.0103±0.0029 

0.0132±0.0026 

0.0118±0.0023 

Tikhonov 

0% 

0.0102±0.0005 

0.0117±0.0003 

0.0117±0.0002 

5% 

0.0102±0.0004 

0.0114±0.0002 

0.0112±0.0001 

10% 

0.0102±0.0003 

0.0108±0.0009 

0.0107±0.0005 

GLS-AC 

0% 

0.01  ±0.001 

0.015±0.0011 

0.0112±0.0003 

5% 

0.0101±0.0014 

0.0146±0.0012 

0.01 06  ±0.0004 

10% 

0.0101±0.0013 

0.0136±0.0009 

0.0112±0.0008 

GLS-LL 

0% 

0.01  ±0.001 

0.0152±0.0012 

0.0113±0.0003 

5% 

0.0101±0.0016 

0.0149±0.0015 

0.0108±0.0006 

10% 

0.0101±0.0016 

0.0138±0.0009 

0.0112±0.001 

Laplacian 

0% 

0.0098±0.0001 

0.0212±0.0001 

0.0112±0.0001 

5% 

0.0098  ±0.0002 

0.0247±0.0001 

0.0097  ±0.0001 

10% 

0.0095  ±0.0001 

0.0276  ±0.0002 

0.0157±0.0128 

Helmholtz 

0% 

0.0099±0.0001 

0.019±0.0002 

0.01 11  ±0.0001 

5% 

0.0099  ±0.0002 

0.0193  ±0.0002 

0.0099  ±0.0001 

10% 

0.0098  ±0.0002 

0.0174±0.0002 

0.0136±0.0001 

Hard-Priors 

0% 

0.0099 

0.0218 

0.0116 

5% 

0.0098 

0.0218 

0.0131 

10% 

0.0098 

0.018 

0.0166 

Table  1(a) 

Method 

Noise 

Region-0 

Region- 1 

Region-2 

level 

Actual 

— 

1.0 

1.0 

3.0 

LM 

0% 

1.0356±0.2364 

0.9995  ±0.0359 

2.3758±0.5160 

5% 

L 075  ±0.0357 

1.0555±0.3254 

1.8215±0.3144 

10% 

1.2672±0.9086 

1.3111  ±0.4128 

1.7111±0.6112 

Tikhonov 

0% 

1.0096  ±0.0397 

1.1153±0.0260 

1.1 644  ±0.0251 

5% 

1.0111  ±0.0004 

1.0912±0.0189 

1.0934±0.0104 

10% 

1.0107±0.0216 

1.0441  ±0.0062 

1.0416±0.0035 

GLS-AC 

0% 

1.0034±0.0688 

1.0335±0.0199 

1.6838±0.1961 

5% 

1.0008±0.0916 

1.0670±0.0362 

1.6972±0.2037 

10% 

0.9987±0.0831 

1.0761  ±0.0343 

1.3703  ±0.0773 

GLS-LL 

0% 

1.0022  ±0.0693 

1.03±0.0183 

1.7801±0.2573 

5% 

0.9998±0.1035 

1.0567  ±0.0329 

1.8502±0.3034 

10% 

0.9981  ±0.0947 

1.0839±0.0425 

1.4271  ±0.0990 

Laplacian 

0% 

0.9918±0.0155 

0.9429±0.0015 

2.8207±0.0491 

5% 

0.9895  ±0.0202 

0.8559±0.0036 

3.693 1  ±  0. 155 1 

10% 

1.0103±0.0124 

0.7447±0.0011 

1.9884±0.0096 

Helmholtz 

0% 

0.9878±0.0154 

1.0518±0.0018 

2.7833±0.0854 

5% 

0.9813±0.0199 

1.1204±0.0081 

3.4252±0.1947 

10% 

0.9884±0.0121 

1.2766±0.01 

2.1761±0.0382 

Hard-Priors 

0% 

0.9919 

0.9266 

2.7332 

5% 

0.9874 

1.0358 

2.345 

10% 

0.9854 

1.3899 

1.822 

Table  1(b) 

=  XLTL.  Moreover,  if  one  assumes  that  A jul=  jm- which  is 
equivalent  to  taking  a  single  step  in  the  iterative  procedure, 
then  Eq.  (19)  maps  into  Eq.  (11)  with  a-  2X.  Hence,  the  LM 
technique  can  be  viewed  as  a  special  case  of  the  Tikhonov 
method,  which  itself  is  a  special  case  of  the  GLS  approach.  It 


is  important,  however,  to  differentiate  LM  from  the  single- 
step  Tikhonov  approach  because  LM  requires  a  to  reach  zero 
asymptotically  with  number  of  iterations,  whereas  in  the 
Tikhonov  scheme,  X  is  constant.  Moreover,  LM  does  not 
involve  parameters  in  the  objective  function. 
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Fig.  4.  Reconstruction  results  (top  of  the  first  row,  abbreviations  are  given 
in  Appendix  A)  are  shown  using  5%  noisy  data  (a)  without  spatial  priors  and 
(b)  with  spatial  priors.  The  top  row  gives  images  of  jma  and  bottom  row 
shows  iis  images. 


Fig.  5.  Reconstruction  results  (top  of  the  first  row,  abbreviations  are  given 
in  Appendix  A)  are  shown  using  10%  noisy  data  (a)  without  spatial  priors 
and  (b)  with  spatial  priors.  The  top  row  gives  images  of  ixa  and  bottom  row 
shows  yu,'  images. 


5.  Stopping  criterion 

The  importance  of  the  stopping  criterion  in  an  iterative 
procedure  cannot  be  ignored.  The  stopping  criterion  used  in 
this  work  is  based  on  the  first-order  conditions  and  data- 
model  misfit,  which  in  the  limit  ensures  that  the  problem  has 
reached  the  global  minima.  The  iterative  procedure  is 
stopped  when  the  L2  norm  of  the  data-model  misfit  (<5)  does 
not  improve  by  more  than  10_10%  or  the  L2  norm  of  the  first 
order  conditions  is  less  than  10_17%.  Beyond  these  values, 
the  round-off  error  dominates.  This  stopping  criterion  is 
more  robust  because  it  involves  first-order  conditions  as 
well. 

IV.  TEST  PROBLEM 

This  section  provides  the  details  of  the  test  problem  con¬ 
sidered  here.  The  optical  property  distributions  used  for  the 
synthetic  data  (y,  noise  added)  generation  are  shown  in  Fig. 
2.  The  diameter  of  the  domain  was  86  mm.  The  background 
optical  properties  were  jma= 0.01  mm-1  and  jul's  =  1.0  mm-1. 
There  were  two  irregular  shaped  targets,  one  in  /jia  with  a 
contrast  of  2:1  to  background  and  one  in  jul's  with  a  contrast 
of  3:1  relative  to  the  background.  A  mesh  consisting  of  4617 
nodes  (corresponding  to  9040  linear  triangular  elements)  was 
used  for  the  generation  of  data.  Sixteen  light  collection/ 
delivery  fibers  were  arranged  equally  spaced  on  the  boundary 
of  the  circle,  where  one  fiber  was  used  as  the  source  while  all 


other  fibers  served  as  detectors  in  turn  which  produced  a  total 
of  240  measurements  [that  is  240  ln(A)  data  points  and  240  0 
data  points].  The  source  was  modeled  as  a  Gaussian  profile 
with  a  full  width  at  half  maximum  of  3  mm  to  represent  the 
light  applied57  and  was  placed  at  a  depth  of  one  transport 
scattering  distance  from  the  tissue  boundary.58  Noise  levels 
of  1%,  3%,  5%,  and  10%  were  added  to  the  modeled  data 
{[ln(A);0]}  to  form  the  experimental  data  (y).  At  the  same 
time,  the  variances  in  the  data  were  also  computed  to  be  used 
in  the  reconstruction  algorithms. 

The  actual  reconstructions  and  forward  modeled  data 
computation  were  performed  on  different  FEM  meshes.59 
This  mesh  has  the  same  diameter  (86  mm)  with  1785  FEM 
nodes,  which  corresponded  to  3418  linear  triangular 
elements.  The  expected  distribution  of  optical  properties  is 
given  in  Fig.  3(a)  (first  column).  Background  optical  proper¬ 
ties  were  used  as  initial  estimates  (jul0)  in  the  evaluation  of 
reconstruction  methods.  The  number  of  parameters  to  be  es¬ 
timated  was  3570  (1785  in  and  1785  in  The  number 
of  data  points  available  for  reconstruction  was  480  [240  of 
ln(A)  and  240  of  0\  The  dimension  of  J  was  480  X  3570,  Ws 
was  480  X  480,  and  was  3570  X  3570.  Optical  prop¬ 

erty  distributions  were  reconstructed  from  the  data  without 
noise  (bias  calculations)  as  well  as  with  noise  levels  of  1%, 
3%,  5%,  and  10%.  The  reconstructions  are  repeated  for  the 
case  of  3%  noise  in  the  data  with  increasing  complexity  (tar¬ 
gets)  in  the  optical  property  distributions. 
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Fig.  6.  A  plot  of  the  rms  error  in  the  estimated  optical  properties  is  shown  as 
a  function  of  increasing  noise  level  for  all  reconstruction  techniques. 


V.  RESULTS  AND  DISCUSSION 

Initially  all  reconstruction  techniques  were  executed  on  a 
data  set  without  noise  to  estimate  the  bias.  Note  that  for  these 
calculations  the  variance  was  found  between  the  data  gener¬ 
ated  using  meshes  (described  in  Sec.  IV)  with  4617  nodes 
and  1785  nodes.  The  results  without  employing  spatial  prior 
information  from  the  reconstruction  techniques  are  given  in 
Fig.  3(a).  The  first  column  shows  the  expected  distribution 
for  the  1785  node  mesh  used  in  the  reconstruction  and  for¬ 
ward  model  calculations.  The  Tikhonov  approach  failed  to 
recover  the  contrast.  This  was  primarily  due  to  the  choice  of 
X,  which  was  based  on  the  maximum  variance  value,  which 
biases  the  problem  to  data  points  that  are  above  the  average 
noise  level.  Since  DOT  is  known  to  have  a  large  dynamic 
range  in  the  data  (at  least  eight  orders  of  magnitude57),  this 
choice  of  X  deemphasize  the  data  points  that  have  low  or 
intermediate  variance  values.  The  root-mean- square  (rms)  er¬ 
rors  between  the  expected  and  reconstructed  optical  proper¬ 
ties  are  plotted  in  Fig.  6.  The  mean  and  standard  deviation  in 
the  reconstructed  images  for  different  regions  [labeled  in  first 
column  of  Fig.  3(a)]  using  the  reconstruction  techniques  dis¬ 
cussed  until  now  are  given  in  Table  I.  In  the  case  of  no 
spatial  priors,  LM  gives  less  bias  in  where  as  GLS  per¬ 
forms  better  in  jul's .  The  bias  calculations  were  repeated  with 
spatial-priors  and  the  reconstruction  results  are  presented  in 
Fig.  3(b).  These  rms  errors  in  the  optical  properties  are  also 
plotted  in  Fig.  6.  Surprisingly  the  soft-prior  approach  (La- 
placian  and  Helmholtz)  performed  better  than  the  hard-prior 
strategy.  It  can  also  be  observed  from  Fig.  6  and  Table  I  that 
the  usage  of  spatial  priors  reduces  the  bias  by  at  least  a  factor 
of  2. 

Figure  4(a)  shows  reconstruction  results  using  data  with 
5%  noise  in  amplitude  without  employing  spatial  priors. 
Once  again  the  Tikhonov  approach  fails  to  recover  the  con¬ 
trast.  The  LM  results  are  dominated  by  boundary  artifacts. 
Figure  4(b)  presents  the  results  from  the  same  data  set  when 
spatial  priors  were  employed.  Figures  5(a)  and  5(b)  show 


similar  kinds  of  effort  for  the  case  of  data  with  10%  noise. 
Note  that,  for  the  same  choice  of  the  regularization  parameter 
(X),  Tikhonov  minimization  scheme  with  spatial  priors 
yielded  quantitatively  more  accurate  results  compared  to 
without  spatial  priors  case,  indicating  that  the  reconstructed 
image  accuracy  along  with  the  quality  largely  depends  on  the 
prior  information.  The  rms  error  in  the  reconstructed  /za  and 
jm's  images  are  plotted  in  Fig.  6  as  a  function  of  increasing 
noise  level.  The  rms  error  using  the  LM  approach  increases 
with  increasing  noise.  GLS  techniques  perform  very  well 
even  in  the  case  of  10%  noise  [Figs.  5(a)  and  6].  Among  the 
GLS  methods,  usage  of  an  analytical  covariance  form  gives 
better  results  (—13%  less  rms  error)  in  /za  and  the  local 
Laplacian  form  performs  slightly  better  (—3%  less  rms  er¬ 
ror)  in  jmfs .  In  the  case  of  employment  of  spatial-priors,  it  can 
clearly  be  seen  [from  Figs.  4(b),  5(b),  and  6  and  Table  I]  that 
hard-priors  perform  better  in  /jl's  reconstruction  when  the 
noise  level  is  below  10%.  Among  the  soft-prior  results,  for 
/za,  the  rms  error  linearly  increases  with  increasing  noise 
level  in  the  Laplacian  case  (Fig.  6).  In  jn's  reconstructions,  the 
performance  of  Laplacian  and  Helmholtz  are  comparable, 
clearly  Helmholtz  performs  slightly  better  (—5%)  when  the 
noise  level  is  above  3%.  Interestingly,  the  Helmholtz  regu¬ 
larization  emerges  with  the  lowest  rms  error  in  fia  recon¬ 
struction.  This  is  primarily  because  of  the  covariance  length 
factor  in  the  Helmholtz  form  of  the  regularization  matrix 
[Eq.  (31)],  which  ensures  that  the  optical  properties  covary 
within  that  correlation  length  (in  here  it  is  5  mm).  The  same 
explanation  is  true  for  the  GLS-analytical  covariance  form 
[Eq.  (28)],  which  performs  better  in  /xa  estimation.  It  is  also 
important  to  note  that  in  the  case  of  a  limited  number  of 
wavelengths,  Srinivasan  et  al .60  have  shown  that  5%  error  in 
the  optical  property  estimate  (juLa  and  jul's)  can  lead  to  ap¬ 
proximately  45%  error  in  spectral  properties  (hemoglobin, 
water  concentrations,  oxygen  saturation,  and  scattering  esti¬ 
mates)  of  tissue.  Any  small  improvement  in  the  optical  prop¬ 
erty  estimates  would  be  important  for  improvement  in  the 
utility  of  this  type  of  imaging  under  practical  conditions. 

To  emphasize  the  effects  of  complexity  on  the  reconstruc¬ 
tion  procedures,  a  set  of  simulations  were  performed  with  an 
increasing  number  of  targets.  Each  target  was  chosen  to  be 
circular  with  a  diameter  of  10  mm.  The  contrast  to  back¬ 
ground  optical  properties  was  2:1.  The  target  locations  and 
corresponding  optical  properties  are  shown  in  the  first  col¬ 
umn  of  Fig.  7(a).  The  targets  were  also  labeled  from  1  to  4 
(background  is  labeled  as  0).  The  data  used  in  this  case  have 
a  noise  level  of  3%.  A  total  of  four  different  reconstructions 
were  performed  by  adding  each  target  at  a  time  (from  1  to  4). 
The  result  of  the  four  target  case  is  shown  in  Fig.  7.  Corre¬ 
sponding  mean  and  standard  deviation  of  the  reconstructed 
optical  properties  for  different  regions  [labeled  in  first  col¬ 
umn  of  Fig.  7(a)]  are  given  in  Table  II.  Figure  8  contains  a 
plot  of  rms  error  in  the  reconstructed  optical  properties  with 
increasing  number  of  targets.  The  rms  error  increases  with 
increasing  number  of  targets  for  every  reconstruction  algo¬ 
rithm.  Note  that  targets  3  and  4  were  placed  close  to  the 
center  of  the  domain,  where  the  sensitivity  is  low  compared 
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Fig.  7.  Reconstruction  results  (top  of 
the  first  row,  abbreviations  are  given  in 
Appendix  A)  are  shown  using  3% 
noisy  data  (a)  without  spatial  priors 
and  (b)  with  spatial  priors  for  four  tar¬ 
gets  in  the  tissue  as  shown.  The  top 
row  gives  images  of  [ia  and  bottom 
row  shows  jul's  images.  The  actual  yufl 
and  jul's  with  target  numbers  are  given 
in  the  first  column  of  (a). 


to  the  periphery.58  Moreover,  increasing  the  /Jia  targets  (from 
1  to  2,  target  numbers  1  and  3),  caused  the  rms  error  to 
increase  by  at  least  30%.  The  same  is  true  with  the  jul's  tar¬ 
gets.  In  the  case  of  multiple  targets,  the  Helmholtz  type  of 
regularization  matrix  resulted  in  the  least  error  in  both 
and  i±s .  Even  though  the  hard-prior  case  performs  very  well 
in  terms  of  lowest  rms  error  for  a  single  target,  as  the  com¬ 
plexity  (or  number  of  parameters  to  be  estimated)  of  the 
problem  increases,  it  clearly  performs  poorer  than  most  of 
the  techniques  presented  (Fig.  8). 

Even  though  the  choice  of  the  Tikhonov  regularization 
parameter  (X)  given  by  Eq.  (14)  is  the  optimal,  the  other 
common  way  is  to  use  L-curve  analysis.61  The  L  curve  for 
DOT  is  shallow,62  similar  to  the  estimation  problem  in  elec¬ 
trical  impedance  tomography,  which  poses  a  problem  in  se¬ 
lection  of  X,  and  is  shown  to  be  unreliable  in  Ref.  59. 

Table  III  gives  the  computational  time  per  iteration  for 
each  of  the  reconstruction  technique  (in  these  two- 
dimensional  cases)  on  Pentium  IV  (dual  core)  2.8  GHz,  2GB 


RAM  Linux  work  station.  GLS  schemes  take  little  more 
computation  time  than  the  Tikhonov  minimization,  as  ex¬ 
pected  hard  priors  took  the  least  computation  time. 

Overall,  the  inclusion  of  spatial  priors  has  an  important 
positive  effect.  The  errors  in  the  estimated  optical  properties 
are  also  reduced  by  at  least  a  factor  of  2  with  spatial  infor¬ 
mation.  The  reconstructed  images  also  contain  the  fine  fea¬ 
tures  extracted  from  conventional  imaging  modalities. 
Through  the  incorporation  of  the  individual  variability  of  the 
data  points  and  optical  parameters  (GLS  scheme),  recon¬ 
struction  performs  better  even  when  the  noise  level  in  the 
data  is  high.  It  is  also  important  to  note  that,  as  mentioned 
before,  iteration  number  8  (which  is  the  best  result  in  terms 
of  lowest  rms  error)  is  chosen  for  rms  error  calculations  in 
LM  approach,  after  this  iteration,  the  solution  becomes  un¬ 
stable.  Whereas  the  rest  of  the  approaches  yield  stable  solu¬ 
tions  (error  in  optical  properties  did  not  increase  with  in¬ 
creasing  iterations).  When  the  individual  data  point  variances 
were  not  considered  (Tikhonov  approach),  the  reconstruction 
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Table  II.  Mean  and  standard  deviation  of  the  reconstructed:  (a)  jma  and  (b)  /x'  values  (in  mm  ')  for  different 
regions  [labeled  in  first  column  of  Fig.  7(a)]  recovered  with  data  having  3%  noise  for  images  shown  in  Fig.  7. 


Method 

Region-0 

Region- 1 

Region-2 

Region- 3 

Region-4 

Actual 

0.01 

0.02 

0.01 

0.02 

0.01 

LM 

0.0101  ±0.0004 

0.0113±0.0001 

0.0112±0.0002 

0.01 11  ±0.0003 

0.011  ±0.0002 

Tikhonov 

0.0102±0.0004 

0.011  ±0.0001 

0.0112±0.0001 

0.0109±0.0001 

0.011  ±0.0001 

GLS-AC 

0.0102±0.0009 

0.0129±0.0003 

0.01 11  ±0.0003 

0.0114±0.0003 

0.0113±0.0003 

GLS-LL 

0.0102±0.0011 

0.0133±0.0004 

0.0115±0.0004 

0.0113±0.0003 

0.0113±0.0002 

Laplacian 

0.01  ±0.0002 

0.0181±0.0001 

0.0105±0.0001 

0.0152±0.0001 

0.0158±0.0001 

Helmholtz 

0.01  ±0.0002 

0.0169±0.0001 

0.0115±0.0001 

0.0149±0.0001 

0.0158±0.0001 

Hard-Priors 

0.01 

0.0158 

0.0126 

(a) 

Region-2 

0.014 

0.0158 

Method 

Region-0 

Region- 1 

Region-3 

Region-4 

Actual 

1.0 

1.0 

2.0 

1.0 

2.0 

LM 

1.0063  ±0.0986 

1.1333±0.0027 

1.24  ±0.0623 

1 . 1 191  ±  0.0396 

1.097  ±0.0366 

Tikhonov 

1.0051  ±0.0217 

1.0341  ±0.0019 

1.0575  ±0.0073 

1.0321  ±0.0056 

1.0329  ±0.0043 

GLS-AC 

0.9993  ±0.0489 

0.9885±0.0139 

1.2486  ±0.0447 

1.021  ±0.0234 

1.1 184  ±0.0076 

GLS-LL 

0.9987±0.0553 

0.9764±0.0127 

1.2726±0.0596 

1.0271  ±0.0262 

1.1422±0.0105 

Laplacian 

0.9886±0.0163 

1.0891  ±0.0023 

2.3799±0.0242 

1.3445  ±0.0043 

1.4044  ±0.0036 

Helmholtz 

0.9899±0.0164 

1.1499±0.0037 

2.1122±0.0386 

1.3382±0.0079 

1.3521  ±0.0066 

Hard-Priors 

0.9856 

1.3712 

1.7319 

(b) 

1.4471 

1.5255 

algorithm  may  not  have  the  ability  to  recover  the  contrast  in 
the  target.  Moreover,  simultaneous  estimation  of  both  ab¬ 
sorption  and  scattering  coefficients  causes  crosstalk  between 
the  two  parameter  estimates.  Even  with  error-free  spatial  pri¬ 
ors,  as  the  complexity  of  the  estimation  problem  (or  number 
of  targets)  increased  for  a  given  noise  level  in  the  data,  the 
parameter-reduction  (hard-priors)  technique  failed  to  give  the 
best  estimates  of  the  optical  properties  due  to  its  LM  imple¬ 
mentation. 


VI.  CONCLUSIONS 

The  diffuse  optical  tomography  inverse  problem  is  often 
solved  by  Levenberg-Marquardt/modified  Tikhonov  minimi¬ 
zation.  A  generalized  approach  for  diffuse  optical  tomogra¬ 
phic  imaging  which  incorporates  the  expected  variability  of 
the  data  noise  and  magnitude  of  the  optical  parameter  varia¬ 
tion  is  presented  as  a  structured  weight-matrix  regularization. 
It  is  also  shown  that  Tikhonov  minimization  and  the 
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2  3 
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Fig.  8.  Plot  of  the  rms  error  in  the  es¬ 
timated  optical  properties  is  shown  for 
increasing  number  of  targets  with  3% 
noise  in  the  data  for  all  reconstruction 
techniques  (legend  of  the  figure).  Ab¬ 
breviations  used  for  the  techniques  are 
given  in  Appendix  A.  The  targets  used 
are  numbered  in  the  images  presented 
in  Fig.  7(a). 
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Table  III.  Comparison  of  computation  time  per  iteration  for  different  recon¬ 
struction  techniques  on  Pentium  IV  (dual  core)  2.8  GHz,  2  GB  RAM  Linux 
work  station,  the  abbreviations  used  for  the  reconstruction  techniques  are 
given  in  Appendix  A. 


Reconstruction  method 

Computation  time  per  iteration 

LM 

17.92  Sec 

Tikhonov 

21.28  Sec 

GLS-AC 

23.39  Sec 

GLS-LL 

23.39  Sec 

Laplacian 

22.78  Sec 

Helmholtz 

22.78  Sec 

Hard-Priors 

10.73  Sec 

Levenberg-Marquardt  approach  are  special  cases  of  this 
GLS  minimization  formalism.  Weight  matrices  that  are  used 
in  this  reconstruction  procedure,  consisting  of  the  variance 
and  covariance  among  the  data  points  and  optical  properties, 
penalize  the  solution  to  match  the  modeled  data  with  the 
experimental  data  more  appropriately.  This  framework  can 
also  be  used  to  incorporate  structural  information,  given  by 
MR,  computed  tomography,  or  other  imaging  modalities 
when  the  two  are  acquired  on  the  same  tissue  volume.  Using 
a  test  problem,  all  of  these  techniques  are  studied  in  terms  of 
the  data  noise  level  and  test  field  complexity  and  a  uniform 
comparison  was  made  using  the  same  implementation 
scheme  for  each  minimization  method.  Even  with  highly 
noisy  data,  the  GLS  approach  gives  meaningful  reconstruc¬ 
tion  results.  It  appears  that  the  standard  Levenberg- 
Marquardt  approach  may  be  unstable  for  the  DOT  problem. 
It  is  also  shown  that  consideration  of  the  individual  variances 
of  data  points  is  the  key  for  an  estimation  procedure  to  re¬ 
cover  high  optical  contrast.  Employing  spatial  information 
reduced  the  errors  in  the  reconstruction  results  by  at  least  a 
factor  of  2.  Parameter  reduction  using  spatial  priors  can  pro¬ 
duce  erroneous  results  when  the  noise  level  is  high.  The 
same  is  true  for  increasing  numbers  of  targets.  Luture  work 
includes  investigating  various  approaches  for  incorporating 
spatial  priors  into  the  GLS  scheme  with  experimental  data 
sets.  Moreover,  a  thorough  examination  of  these  techniques 
in  three-dimensional  case  will  be  taken  up  as  a  future  inves¬ 
tigation.  The  computer  algorithms  and  test  data  used  in  this 
article  (along  with  some  additional  material)  are  given  at  this 
web  page.63 
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APPENDIX  A:  TERMINOLOGY 

DOT — diffuse  optical  tomography. 

jma(r) — optical  absorption  coefficient  of  tissue. 


jul's  (r) — reduced  (or  transport)  scattering  coefficient  of  tis¬ 
sue. 

D( r) — optical  diffusion  coefficient  of  tissue=  l/3[/xa(r) 
+  /^(r)]. 

fJL — parameters  (generalized)  to  estimate =[D(r) 
julq — prior  value  of  the  parameters  (initial  guess,  generally 
obtained  from  prior  calibration  of  data45,46). 

F(hjl) — forward  model. 

G(/ul) — Modeled  data  (G — sampled  forward  model 
=S{F}). 

A — amplitude  of  the  signal. 

8 — phase  of  the  signal. 
y — Measured  data=[ln(A) ;  6\ 

||X|| — L2  norm  of  vector  X=^f=1X2. 

& — data-model  misfit =y-G(jm). 

Ws — weight  matrix  for  £=[cov(£)]-1  (Appendix  A-4  of 
Ref.  47). 

weiSht  matrix  for  ya-/x0=[cov(/x-/x0)]“1  (Ap¬ 
pendix  A-4  of  Ref.  47). 

\ — Tikhonov  regularization  parameter. 

L — Tikhonov  regularization  matrix. 

I — identity  matrix. 
o2 — variance. 

J — Jacobian  of  the  sampled  forward  model =dG(jm)/  djm. 
fl — objective  function. 

Error — true  value-estimated  value  (prediction). 

Bias — difference  between  the  true  optical  property  distri¬ 
bution  and  estimated  optical  properties  in  the  case  of  model 
generated  data  (without  adding  the  noise). 

Ill-posed — Small  changes  in  the  data  can  cause  large 
changes  in  the  parameters. 

Ill-conditioned — the  condition  number  (ratio  of  largest 
singular  value  to  smallest  singular  value)  is  large,  which  im¬ 
plies  the  inverse  solution  would  not  be  unique. 

Ill-determined — (or  under-determined)  the  number  of  in¬ 
dependent  equations  are  smaller  than  number  of  unknowns. 
Unstability — error  gets  amplified  with  iterations. 

LM — Levenberg-Marquardt  minimization  (Sec.  Ill  A  1). 
Tikhonov — Tikhonov  minimization  scheme  without 
spatial-priors,  L-I  (Sec.  Ill  A  2). 

GLS-AC — generalized  least  squares  minimization 

scheme  (Sec.  Ill  A  3)  with  analytical  covariance  form  for 
[Eq.  (28)]. 

GLS-LL — generalized  least  squares  minimization 

scheme  (Sec.  Ill  A  3)  with  local  Laplacian  form  for 
[Eq.  (29)]. 

Laplacian — Tikhonov  minimization  scheme  in  the  case 
of  soft  priors  (Sec.  Ill  B  1)  where  L  approximates  Laplacian 
form,  defined  by  Eq.  (30). 

Helmholtz — Tikhonov  minimization  scheme  in  the  case 
of  soft  priors  (Sec.  Ill  B  1)  where  L  approximates  Helmholtz 
form,  defined  by  Eq.  (31). 

Hard-priors — parameter-reduction  technique  based  on 
spatial  priors  (Sec.  Ill  B  2). 
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APPENDIX  B:  TIKHONOV  REGULARIZATION- 
SINGULAR  VALUES 

It  is  interesting  to  examine  Tikhonov  regularization  from 
the  point  of  view  of  singular  values.  If  one  rewrites  the  up¬ 
date  equation  [Eq.  (19)]  as 

[JTJ  +  X/]A/*f  =  JT  St_i  +  C,  (B 1) 

where  C=X(/z/_1-/z0),  as  it  is  a  constant  vector  for  a  chosen 
iteration  i.  Singular- value  decomposition  of  J  gives 

J=USVT ,  (B2) 


where  U  and  V  are  orthonormal  matrices  containing  the  sin¬ 
gular  vectors  of  /,  i.e.,  UTU=I  and  VTV=I.  S'  is  a  diagonal 
matrix  containing  the  singular  values  (St)  of  J.  Substituting 
this  into  update  equation  [Eq.  (B 1)]  generates 

C VSTUTUSVT+  XI)  Ajuit  =  VSTUTSi_i  +  C.  (B3) 

Using  the  orthonormal  properties  of  U  and  left  multiply¬ 
ing  by  VT  on  both  sides  of  Eq.  (B3)  yields 

C VTVSTSVT  +  XV7)  A  M  =  VTVSTUTSi_]  +  VTC.  (B4) 

Now  using  the  orthonormal  properties  of  V  and  rearrang¬ 
ing  the  terms  leads  to 

(STS  +  XI)  VTAjmi  =  STUTSi_l  +  VTC.  (B5) 


Taking  the  inverse,  left  multiplying  by  V  and  simplifying 
the  result  gives 

A  &  =  V(STS  +  Xiy'iS^S^  +  VTC) .  (B6) 


Writing  Eq.  (B7)  in  the  form 

A  jULi  =  VDP,  (B7) 


where  P = {STUT8i_  x  +  VTC) ,  a  column  vector,  and  D  is  a  di¬ 
agonal  matrix  which  has  the  form 


0 


^  +  x 


if  i  *  j 


(B8) 


Similar  expressions  hold  for  L¥=I  (Ref.  65)  in  Eq.  (18). 
Considering  the  case  X  =  0,  one  can  clearly  see  that  for  an 
ill-conditioned  matrix  J ,  implying  some  of  the  singular  val¬ 
ues  are  almost  zero  (5)~  0),  the  inversion  becomes  unstable 
(some  of  the  diagonal  values  of  D  become  infinite).  By  using 
Tikhonov  regularization,  even  when  5)  =  0,  the  inversion  pro¬ 
cedure  is  stabilized  [Eq.  (B8)].  The  X  act  as  a  filtering  factor, 
giving  the  name  Tikhonov  filtering65  for  this  procedure. 
Moreover,  as  this  X  damps  the  amplification  of  the  diagonal 
values  of  D  for  smaller  values  of  St  in  Eq.  (B8),  this  is  also 
known  as  damped  least  squares  minimization  procedure.65 
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Abstract:  Near-Infrared  (NIR)  tomographic  image  reconstruction  is  a 

non-linear,  ill-posed  and  ill-conditioned  problem,  and  so  in  this  study,  dif¬ 
ferent  ways  of  penalizing  the  objective  function  with  structural  information 
were  investigated.  A  simple  framework  to  incorporate  structural  priors 
is  presented,  using  simple  weight  matrices  that  have  either  Laplacian  or 
Helmholtz-type  structures.  Using  both  MRI-derived  breast  geometry  and 
phantom  data,  a  systematic  and  quantitative  comparison  was  performed 
with  and  without  spatial  priors.  The  Helmholtz-type  structure  can  be  seen 
as  a  more  generalized  approach  for  incorporating  spatial  priors  into  the 
reconstruction  scheme.  Moreover,  parameter  reduction  (i.e.  hard  prior 
information)  in  the  imaging  field  through  the  enforcement  of  spatially 
explicit  regions  may  lead  to  erroneous  results  with  imperfect  spatial  priors. 
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1.  Introduction 

Near  Infrared  (NIR)  optical  tomography  is  a  method  that  uses  light  in  the  650-900nm  wave¬ 
length  range  to  recover  images  of  the  internal  spatial  distribution  of  tissue  optical  properties, 
absorption  (or  chromophore  concentrations)  and  scattering  parameters  [1-3].  When  imaged  at 
multiple  wavelengths,  these  quantities  can  be  used  to  estimate  tissue  hemoglobin  and  water 
concentration^]  and  have  the  advantage  of  being  acquired  non-invasively  and  without  ionizing 
radiation.  The  imaging  procedure  can  be  rapidly  or  repeatedly  applied  to  investigate  physio¬ 
logical  state,  and  systems  can  be  integrated  into  conventional  imaging  platforms  such  as  X-ray 
mammography,  Ultrasound  and  MRI[4-9].  These  hybrid  systems  have  been  shown  to  achieve 
superior  performance  in  terms  of  resolution  and  quantitative  accuracy  which  should  provide 
more  accurate  physiological  data  from  the  tissue  under  investigation^- 16].  However,  a  funda¬ 
mental  question  is  how  to  utilize  the  spatial  information  from  the  clinical  system  to  maximize 
the  accuracy  of  NIR  tomography.  In  this  study,  the  ability  to  improve  the  quantitative  accuracy 
of  regions  imaged  with  NIR  tomography  was  investigated,  in  the  setting  where  prior  spatial 
information  is  readily  available. 

This  work  explores  image  reconstruction  strategies  that  take  advantage  of  multimodality  im¬ 
age  data,  in  particular,  the  combination  of  MRI  with  NIR  optical  tomography  for  breast  cancer 
imaging.  MRI  provides  structural  information  at  high  spatial  resolution  (near  1  mm),  whereas 
NIR  imaging  has  relatively  poor  spatial  resolution  (near  4-7  mm).  Yet  MR  imaging  would  ben¬ 
efit  from  the  molecular- specific  signatures  available  through  NIR[3,  16-18],  specifically  tissue 
hemoglobin  content,  oxygenation  level,  and  water,  as  well  as  scattering  particle  size  and  number 
density  [19].  The  inverse  problem  (image  reconstruction  procedure)  in  NIR  imaging  is  known 
to  be  a  non-linear,  ill-posed  and  ill-conditioned  problem  [20].  Use  of  structural  information  in 
NIR  reconstruction  schemes  has  been  explored  by  several  research  groups.  For  example,  Li  et. 
al  [7]  have  used  the  data  derived  from  X-ray  mammography  for  choosing  different  regulariza¬ 
tion  parameters  for  the  region  of  interest  (ROI)  and  surrounding  tissue,  and  have  shown  that 
the  contrast  and  resolution  of  the  reconstructed  images  can  be  improved.  Srinivasan  et  al.  [21] 
have  developed  a  three-step  reconstruction  process  for  improving  the  quantification  accuracy  of 
small-objects  in  NIR  tomography,  where  they  use  the  conventional  NIR  reconstructed  images 
(first  step)  as  a  structural  prior  for  the  last  two  steps.  Earlier  papers  have  shown  that  optical  con¬ 
trast  can  be  correlated  to  MR  contrast  [6,  9,  13]  and  that  structural  MRI  images  can  be  used  to 
reduce  the  number  of  unknown  parameters  to  be  estimated  [14].  The  difficulty  with  parameter 
reduction  approaches  (referred  to  as  hard-priors  [22])  is  the  potential  of  introducing  error  by 
imposing  incorrect  model  assumptions  and  introducing  variation  due  to  uncertainty  in  the  prior 
information  (even  when  the  underlying  model  is  appropriate).  For  example,  the  features  which 
lead  to  contrast  in  one  imaging  system  may  not  be  spatially  coregistered  with  those  that  pro¬ 
duce  contrast  in  another  imaging  system.  Further,  segmentation  of  congruent  features  always 
includes  classification  errors  resulting  from  digitization.  Recently,  Boverman  et.  al  [10]  showed 
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that  even  imperfect  priors  which  encode  breast  background  structure  improves  anomaly  local¬ 
ization,  but  at  the  expense  of  biasing  the  spectroscopic  dimension  of  the  image  reconstruction. 
Another  type  of  approach  for  constraining  the  problem,  often  described  as  soft-priors,  penal¬ 
izes  the  variation  within  regions  which  are  assumed  to  have  the  same  properties  by  controlling 
regularization.  Brooksby  et.  al  [15,  16,  23]  have  developed  a  Laplacian  type  of  regularization 
that  allows  intra-region  variability.  This  is  a  method  which  works  well  even  if  the  confidence  in 
the  prior  structural  information  is  low. 

This  paper  develops  a  more  generalized  framework,  known  as  soft-priors,  for  incorporating 
the  structural  priors  into  the  NIR  image  reconstruction  process,  and  explores  a  covariance-based 
constraint  scheme  adopted  from  finite  differencing  of  the  Helmholtz  equation  in  addition  to  the 
soft  and  hard  prior  approaches  noted  above.  The  soft-prior  approach  allows  optical  property 
variation  with  a  given  region,  reducing  biases  caused  by  the  use  of  imperfect  prior  information. 
The  results  indicate  that  imperfect  structural  information  can  generate  errors  in  the  hard-priors 
case,  whereas  the  soft-priors  are  able  to  quantify  regions  more  appropriately.  Simulation  and 
experimental  studies  are  performed  to  demonstrate  the  superior  reconstruction  image  quality. 
These  types  of  procedures  are  needed  to  improve  NIR  imaging  both  in  terms  of  high  spatial 
resolution  available  from  MRI  and  high  contrast  inherent  in  the  NIR  signal. 

2.  Methods:  mathematical  formulation 

2.7.  Diffusion-based  Light  Transport  Model 

Light  transport  in  breast  tissue  can  be  described  accurately  by  the  Diffusion  Equation  (DE), 
which  is  an  approximation  to  the  Radiative  Transport  Equation  (RTE)  [24].  In  the  frequency- 
domain,  the  DE  is  given  by 

— V.D(r)V<t>(r,  (o)  +  (jUa(r)  +  ico/c)<$>(  r,  co)  =  Q0(  r,  co)  (1) 

where  <f>(r,  co)  is  the  photon  density  at  position  r  and  light  modulation  frequency  is  given  by  co 
(in  this  work,  co  =  100  MHz).  The  isotropic  source  term  is  represented  by  Q0( r,  co)  and  speed  of 
light  in  tissue  by  c.  jUa(r)  is  the  optical  absorption  coefficient  and  D(r)  is  the  optical  diffusion 
coefficient,  which  is  defined  as 


D(r)  3[Ma(r)+At'(r)]  (2) 

where  /i'(r)  is  the  reduced  scattering  coefficient,  which  is  defined  as  /l'  =  /is(l  —g).  ps  is  the 
scattering  coefficient  and  g  is  the  anisotropy  factor.  A  Robin-type  (Type-III)  boundary  condition 
is  applied  to  exactly  model  the  refractive-index  mismatch  at  the  boundary  [25].  The  boundary 
data  for  a  frequency  domain  system  are  the  amplitude  and  phase  of  the  measured  signal,  which 
is  used  with  a  Finite  Element  Method  (FEM)  based  reconstruction  procedure  to  obtain  the 
internal  spatial  distributions  of  pa  and  fl's. 

2.2.  Standard  image  reconstruction 

The  objective  function  (f 1 )  for  this  procedure  can  be  written  as 

^  Dm“a  {||y-F(D,Aia)||2  +  A||(D,Ata)-(D0,Aia0)||2}  (?) 

Where,  F  is  the  forward  operator  that  generates  the  model  response  and  y  is  the  experimental 
measured  data.  ||.||2  represents  L2-Norm  of  the  vector.  This  is  also  known  as  the  Tikhonov 
approach  [27],  where  A  is  the  regularization  parameter  that  balances  the  current  estimate  of 
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optical  properties  and  the  initial  values,  and  the  data-model  misfit.  More  specifically,  it  is  the 
ratio  of  the  variances  of  the  data  noise  and  the  parameter  (A  =  Oy/o^D  ^)[28].  Minimizing 
Eq.  (3)  (by  setting  first  derivatives  with  respect  to  {ia  and  D  to  zero)  leads  to  an  update  equation 

(. 'TJ  +  A/) (5D,  a)  =  JT(y-  F(D,  ;Ua))-A [(D, *ta)  -  (Do, Mao)]  (4) 

Note  that  deriving  this  update  equation  is  not  so  straight  forward,  the  full  derivation  is  given 
in  Ref.  [28].  In  Eq.  4,  J  is  the  Jacobian  matrix  and  I  is  the  identity  matrix.  Note  that  JTJ  is 
ill-conditioned;  A/  stabilizes  the  matrix.  However,  a  slight  deviation  from  this  update  equation 
is  generally  employed,  which  is  also  known  as  the  Levenberg-Marquardt  (LM)  regularization 
method[29,  30],  assuming  [(<5Z),<5/4)  =  (D,/ia)  —  (Do,/iao)]  leading  to  [31] 

(JT  J  +  2A/)  (<5D,  <5jUa)  =JT(  y-F(D,Aia))  (5) 

Most  of  the  literature  reports  A*  =  2 A  [20,  31],  which  is  true  only  for  the  Levenberg-Marquardt 
minimization  which  does  not  involve  the  parameter  field  in  the  objective  function  (Eq.  3) [29, 
30].  A  detailed  discussion  about  different  least-squares  minimization  methods  is  presented  in 
Ref.  [28].  In  this  LM  approach,  A*  typically  starts  being  the  ratio  of  the  variances  and  is  reduced 
at  each  of  the  iterations  by  a  small  factor  (in  here,  it  is  \/T0  and  also  multiplied  by  the  maximum 
of  the  diagonal  values  of  JTJ[ 21]).  The  iterative  procedure  is  repeated  until  experimental  data 
matches  with  modeled  data  within  a  preset  value  £  (^  data  noise  level).  In  general,  the  initial 
values,  (Dq,  Jtao),  are  obtained  from  a  pre-reconstruction  step  where  the  data  is  calibrated  by  a 
homogeneous  fitting  procedure[32,  33]. 

23.  Inclusion  of  a  priori  information:  Soft-Priors 

The  objective  function  with  inclusion  of  prior  information  is  given  as [15] 


D™  tlly-F(D,Jua)||2  +  A||Z.[(D,Jua)-(D0,JUao)]||2}  (6) 

Here  also  A  is  the  ratio  of  the  variance  of  the  data  noise  to  parameter  field  and  L  is  a  penalty 
matrix  (dimensionless  in  all  the  cases  considered  in  this  work)  which  can  be  derived  from  MRI 
structural  priors  as  indicated  below.  The  update  equation  resulting  from  this  procedure  is: 

(JTJ  +  ALrL)(5D,5/ia)  =  7r(y  —  F(D,/ia))  —  ALrL[(D,^a)  —  (D0,/iao)]  (7) 


In  this  work,  each  location  in  the  computational  discretized  model  is  labeled  according  to  tissue 
type  (fatty,  fibroglandular  or  tumor)  determined  from  MRI  T1 -weighted  images[15,  16,  23].  It 
was  also  assumed  that  there  is  no  covariance  between  the  different  regions  of  the  imaging 
domain.  Since  the  domain  model  does  not  itself  change  throughout  the  iterative  reconstruc¬ 
tion  algorithm,  the  L-matrix  is  calculated  before  the  reconstruction  procedure  and  it  is  used 
through  out  the  process  to  penalize  the  solution.  In  implementation  of  Eq.  7  for  simultaneous 
reconstruction  of  D  and  fj,a,  using  the  Jacobian  form  described  in  Ref.  [18]  leads  to  the  block 
matrices 


HD,a  1 

hi 

r  &d)ltl 

0  1 

H Djla 

H ^  \ 

[  0 

1  (V)i'i  J 
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8k 


(JT)D(y-F(J)^a)) 


(AMjLrL(Ma,_i-jUao)  J 


where  H  represents  Hessian  matrix  ( JTJ ).  Since  D  and  jla  are  considered  to  be  independent 
of  each  other  in  the  estimation  procedure,  the  cross-terms  in  the  combined  regularization  ma¬ 
trix  (A LT L)  become  zero.  Two  forms  for  the  L-matrix  are  discussed  in  the  subsections  below, 
including  the  Laplacian  and  Helmholtz  structures. 
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2.3.1.  a  Laplacian-type  structured  regularization  matrix 
The  Laplace  equation  in  parameter  u(r)  can  be  written  as 

V2w(r)  =  0  (9) 


A  finite  difference  approximation  to  this  equation  for  ‘N’  number  of  equi-space  (step  size  =  h) 
nodes  can  be  written  as  [34] 

V2u(r)h2  «  u\  +  U2  +  . . .  —Nun/ 2  +  •  •  .  +  un-i  +  un  =  0  (10) 

Dividing  the  whole  equation  by  ‘-N’  leads  to 


—u\ 

IsT 


—  U2 

~N~ 


+  •  •  •  +  Un/2  +  •  •  •  + 


—un-  i 
N 


(ID 


The  L  matrix  is  a  matrix  that  relates  each  nodal  property  of  the  numerical  model  to  all  other 
nodes.  Therefore  given  a  node  i  within  the  mesh,  its  relationship  to  another  node  j  having 
Laplacian  structure  (Eq.  11)  within  the  same  mesh  can  be  given  as  [15,  16] 


{0  if  i  and  j  are  not  in  the  same  region 

—  l/N  if  i  and  j  are  in  the  same  region  (12) 

1  if  i=j 


where  N  is  the  number  of  finite  element  mesh  nodes  comprising  a  given  region.  In  this  case, 
LT L  approximates  a  second-order  Laplacian  smoothing  operator  within  each  region,  and  func¬ 
tionally  works  to  average  the  update  within  a  region,  while  allowing  discontinuity  between 
different  regions. 


2.3.2.  a  Helmholtz- type  structured  regularization  matrix 

The  Helmholtz  equation  in  parameter  u(r)  for  a  damped  wave  can  be  written  as 

V2«(r)  -  K2u(r )  =  0  (13) 


where  k  is  the  wave  number,  specifically  k  =  1//,  where  l  covariance  length[34].  I  also  repre¬ 
sents  the  decay  length  scale  over  which  the  parameter  u(r)  has  correlation.  Making  K  =  0,  will 
give  the  Laplace  equation  (Eq.  9).  A  finite  difference  approximation  to  this  equation  for  ‘N’ 
number  of  equi-space  (step  size  =  h)  nodes  can  be  written  as  [34] 

(V2  -  k2)  u(r)h 2  «  u\  +  W2  + . . .  +  [—(N+(ich)2)]uN/2  + . . .  +  un-i  +un  =  0  (14) 

Dividing  the  whole  equation  by  ‘  —  (N  +  ( /c/z)2)  ’  gives 


U\ 

-(N+(Kh)2) 


U2  Un— 1 

-(jV+(K/*)2)  +  ' ' '  +  UN/2  +  ■  ■  ■  +  _  (V+  (ff/i)2) 


Writing  this  as  a  generalized  L-matrix  form  similar  to  Eq.  12 


~{N+{Kh)2) 

(15) 


Uj  = 


0 

1 

N-\-(ich)2 

1 


if  i  and  j  are  not  in  the  same  region 
if  i  and  j  are  in  the  same  region 
ifi=j 


(16) 


For  the  FEM  nodes  case,  h  is  chosen  to  be  the  distance  between  the  nodes.  Moreover,  k  =  l/l 
is  generally  chosen  to  be  the  inverse  of  the  size  of  the  feature  (tumor  in  this  case)  in  the  imaging 
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domain.  This  can  be  seen  as  using  the  best  priors  for  the  estimation  of  optical  properties [28], 
termed  as  best  priors  estimate  (BPE).  As  the  prior  structural  information  is  available  through 
MRI,  /  is  chosen  to  be  the  diameter  of  the  target  (tumor)  in  the  BPE  case.  In  this  case,  LT L  (L 
given  by  Eq.  16)  approximates  a  second-order  Helmholtz  smoothing  operator.  To  determine  the 
effect  of  K  on  the  parameter  reconstruction,  different  values  for  K  are  chosen.  It  is  shown  that 
for  small  values  of  K ,  which  corresponds  to  a  large  correlation  length  (/),  both  Laplacian  and 
Helmholtz  structures  recover  the  same  optical  property  distribution. 


2.4.  Inclusion  of  a  priori  information:  hard-priors 

Reduction  of  parameter  space  to  the  number  of  regions  segmented  from  a  high  resolution  imag¬ 
ing  modality  is  known  as  hard-priors  [2 8].  In  this  section,  a  detailed  formulation  for  this  ap¬ 
proach  is  given.  The  estimation  optical  properties  ( D  and  pa)  in  this  procedure  is  performed 
through  LM  minimization,  which  was  discussed  in  detail  in  Ref.  [28].  The  update  equation 
in  this  case  is  given  by  Eq.  5,  with  Jacobian  having  dimension  of  (2*NM)x(2*NR)  instead  of 
(2*NM)x(2*NN).  In  here,  NM,  NN  and  NR  represent  the  number  of  measurements,  number  of 
FEM  nodes,  and  number  of  regions,  respectively.  The  multiplication  factor  2  in  the  dimensions 
results  from  the  treatment  of  amplitude  and  phase  separately  for  the  frequency-domain  signal, 
and  the  estimated  parameters  being  D  and  pa.  This  transformation  requires  multiplication  of 
original  Jacobian,  /,  by  region  mapper  matrix  (R)[ 22]. 


J  =  JR 


(17) 


where  the  R  has  dimension  of  (2*NN)x(2*NR),  having  a  block  form  R  = 


£ 

£ 


.  J  can  be  also 


seen  as  a  block  matrix  with  J  =  [Jo  /uj.  If  RLj  represented  a  segmented  region  (j)  in  the  FEM 
mesh  for  j  =  1,  2,  ...NR,  then  £  has  the  from 


Zy  = 


if  i  G  RLj 
otherwise 


(18) 


Overall,  the  new  Jacobian  /  matrix  elements  were  produced  by  adding  the  sensitivity  of  nodes 
belonging  to  the  same  region  [22].  Using  /  instead  of  /,  in  Eq.  5  gives  the  update  of  optical 
properties  of  the  segmented  regions  ^<5D,  <5/ia^  having  dimension  of  (2*NR)xl.  To  map  back 
this  solution  vector  to  the  original  dimensions  (2*NN)xl,  requires 

(5D,5M.)  =  *(5D,5i4)  (19) 

It  is  also  important  to  relaize  that  even  though,  NR<CNM,  solving  this  might  still  be  ill- 
posed[20],  leading  to  LM  update  equation[28]. 


3.  Methods:  simulations  and  experiments 

3.1.  Breast  geometry  -  Effect  of  imperfect  a  priori  information 

The  techniques  described  in  Sec.  2  were  used  to  reconstruct  images  from  synthetic  data  with 
1%  random  noise  added.  Numerical  experiments  using  synthetic  data  generated  on  MRI  Tl- 
weighted  breast  images  with  incorrect  priors  to  show  the  effectiveness  of  soft-priors.  Figure 
1(a)  (first  column)  shows  the  original  distribution  of  three  tissue  layers,  namely  fatty  (pa  = 
0.006  mm-1  and  p's  =  0.6  mm-1),  fibro  glandular  ( jLla  =  0.012  mm-1  and  p's  =  1.2  mm-1),  and 
tumor  (pa  =  0.018  mm-1  and  p's  =  1.8  mm-1)  for  the  breast  geometry  (labeled  as  0,  1  and  2 
respectively,  in  Fig.  1(a)  first  column).  Sixteen  light  collection/delivery  fibers  were  arranged 
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Actual  Laplacian  Helmholtz  Hard  Priors 


(a) 


(b) 

Fig.  1.  (a)  Simulated  jia  and  fi's  distributions  from  a  breast  (obtained  from  a  volunteer)  are 
shown  in  the  first  column.  Optical  properties  for  the  region  labeled  ‘O’  (fat)  are:  i^a  =  0.006 
mm-1  and  fis  =  0.6  mm-1.  Region  ‘1’  (fibroglandular)  values  are:  fia  =  0.012  mm-1  and 
fi's  =  1.2  mm-1.  Region  ‘2’  (tumor)  values  are:  /ia  =  0.018  mm-1  and  fd's  =  1.8  mm-1. 
Reconstructed  jia  and  /i'  images  from  different  techniques  with  simulated  data  having  1% 
random  noise  and  imperfect  structural  information  in  defining  region  ‘1’  (7%  reduction 
compared  to  the  original  segmentation)  are  shown  in  the  rest  of  the  columns.  The  middle 
two  columns  use  soft  prior  structural  information  while  the  last  column  shows  the  result 
with  hard  prior  information.  In  the  Helmholtz  case,  K  =  1/8  mm-1  (BPE)  was  used,  (b) 
Cross-sectional  plots,  along  the  dotted  line  in  the  actual  image  (see  first  column  of  (a)),  of 
true  and  reconstructed  jda  and  /r'  distributions. 
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equally  spaced  on  a  circle  (indentions  in  Fig.  1(a)).  In  succession,  one  fiber  was  used  as  the 
source  while  all  other  fibers  served  as  detectors  which  provided  a  total  of  240  measurements.  In 
these  studies,  the  source  was  modeled  as  a  Gaussian  profile  with  a  Full  Width  Half  Maximum 
(FWHM)  of  3  mm  to  most  accurately  represent  the  light  applied  in  the  clinical  system[35],  and 
was  placed  at  a  depth  of  one  transport  scattering  distance  from  the  tissue  boundary.  A  mesh  of 
1785  nodes  (corresponding  to  3418  linear  triangular  elements)  with  uniform  nodal  density  was 
used  for  the  diffusion  model  and  reconstruction  calculations [3 6].  A  total  of  7%  of  the  glandular 
layer  (label- 1)  FEM  nodes  were  labeled  (relative  to  the  original  glandular  layer  nodes)  as  fat 
(label-0)  to  introduce  imperfect  structural  priors. 

The  same  initial  estimates  (optical  properties  of  region  ‘0’)  were  used  as  homogeneous 
starting  conditions.  The  iterative  procedure  was  stopped  once  the  data- model  misfit  (residual) 
did  not  improve  by  more  than  2%  when  compared  with  the  previous  iteration.  The  starting 
value  for  A  is  chosen  to  be  25000  and  75  for  iia  and  D  respectively,  derived  from  the  noise 
characteristics [28],  for  Eq.  5.  For  the  case  of  Eq.  7,  a  starting  value  of  10  was  chosen  and  was 
decreased  as  the  iterations  progressed,  this  procedure  is  outlined  in  Ref.  [21]. 


Table  1.  Mean  and  standard  deviation  of  the  reconstructed  (a)  jj,a  and  (b)  fi's  values  in  differ¬ 
ent  regions  (labeled  in  first  column  of  Fig.  1(a))  recovered  with  simulated  data  having  1% 
random  noise  and  imperfect  structural  information  defining  region  T’  (7%  reduction  com¬ 
pared  to  the  original  segmentation).  The  corresponding  reconstructed  images  are  shown  in 
Fig.  1(a) 


Methods 

Region-0 

Region-1 

Region-2 

Actual 

0.006 

0.012 

0.018 

Laplacian 

0.0064±0.0010 

0.0117±0.0018 

0.0156±0.0018 

Helmholtz  (k  =  1/8) 

0.0062±0.0011 

0.0120±0.0020 

0.0156±0.0017 

Hard  Priors 

0.006 

0.0118 

0.0843 

(a) 


Methods 

Region-0 

Region- 1 

Region-2 

Actual 

0.6 

1.0 

1.8 

Laplacian 

0.63±0.10 

1 . 13±0. 1 8 

1.67±0.23 

Helmholtz  (k  =  1/8) 

0.64±0.09 

1.09±0.16 

1.64±0.22 

Hard  Priors 

0.64 

1.13 

0.23 

(b) 


3.2.  Breast  geometry  -  effect  of  data  noise  level 

The  techniques  described  in  Sec.  2  were  used  to  reconstruct  images  from  synthetic  data  with  1, 
3,  5  and  10%  Gaussian  distributed  noise  to  see  the  effect  of  data  noise  level  on  the  reconstruction 
techniques.  The  breast  geometry  (and  the  optical  properties)  were  equivalent  to  the  previous 
section  (Fig.  1(a)),  however,  perfect  spatial  priors  were  used.  The  same  FEM  mesh  as  described 
above  was  employed  in  the  forward  and  reconstruction  problems.  Optical  properties  of  region 
‘O’  were  used  as  initial  guess  for  the  iterative  procedure.  The  regularization  parameter  (A)  and 
stopping  criterion  was  chosen  according  to  each  data  noise  level. 

3.3.  Phantom  studies 

A  multi-layered  gelatin  phantom  (86  mm  diameter,  25  mm  height)  was  fabricated  with  different 
optical  properties  using  heated  mixtures  of  water  (80%),  gelatin  (20%)  (G2625,  Sigma  Inc.), 
India  ink  for  absorption,  and  Ti02  (titanium  oxide  powder,  Sigma  Inc.)  for  scattering[15]  (see 
Fig.  3).  Different  layers  of  gelatin  were  constructed  by  successively  hardening  gel  solutions 
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Fig.  2.  (a)  Reconstructed  [ia  and  ji's  images  from  different  techniques  with  simulated  data 
having  5%  random  noise  and  perfect  structural  priors  (actual  images  are  shown  in  the  first 
column  of  Fig.  1(a)).  The  first  column  shows  the  reconstruction  results  without  the  use  of 
prior  information.  The  middle  two  columns  use  soft  prior  structural  information  while  the 
last  row  shows  the  result  with  hard  prior  information.  In  the  Helmholtz  case,  K  =  1/8  mm-1 
(BPE)  was  used,  (b)  The  mean  values  and  standard  deviations  (plotted  as  error  bars)  in  jj,a 
and  ii's  for  different  regions  of  breast  geometry  (labeled  in  actual  image)  with  increasing 
noise  level  (1%  to  10  %). 
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containing  different  amounts  of  ink  and  T1O2.  A  cylindrical  hole  (diameter  of  16  mm  and 
height  of  24  mm)  was  filled  with  liquid  to  mimic  the  tumor.  The  first  column  in  Fig.  4  shows 
the  axial  cross-section  of  three-layers  of  the  phantom  (Fig.  3)  where  the  region  labeled  ‘0’  has 
the  optical  properties  (jda  =  0.0065  mm-1  and  fl's  =  0.  65  mm-1),  similar  to  the  typical  fatty 
layer  in  the  breast [6]  and  a  thickness  of  10  mm.  The  fibroglandular  layer  (diameter  76mm) 
has  optical  properties  (region  labeled  ‘1’)  of  ila  =  0.01  mm-1  and  /l'  =  1.0  mm-1.  The  tumor 
(represented  by  the  region  labeled  ‘2’)  has  a  diameter  of  16  mm  with  optical  properties  of 
fj,a  =  0.02  mm-1  and  fl's  =  1.2  mm-1.  The  optical  properties  were  validated  by  measuring 
large  cylindrical  samples  of  each  layer.  Appropriate  mixtures  of  Intralipid  and  India  ink  were 
used  to  achieve  the  desired  optical  parameters  of  the  tumor.  Data  was  acquired  using  a  clinical 
NIR  system [3 5]  where  the  fibers  were  marked  and  photographed  to  extract  region  information 
(analogous  to  MRI  images).  This  regional  information  was  used  to  label  the  corresponding 
regions  in  the  FEM  mesh[15].  A  mesh  of  1785  nodes  (corresponding  to  3418  linear  triangular 
elements)  was  used  for  the  diffusion  model  calculations  and  a  mesh  having  1360  nodes  was 
used  in  the  reconstruction[36].  The  meshes  considered  in  this  work  have  uniform  nodal  density 
across  the  total  computation  volume.  NIR  data  was  calibrated  using  a  reference  homogenous 
phantom  to  obtain  initial  optical  properties  estimates  and  minimize  the  variation  between  the 
16  optical  channels  according  to  standard  practice  in  human  imaging  studies [32,  33]. 


Fig.  3.  Photograph  for  gelatin  phantom  (representing  the  idealized  two-dimensional  cross- 
sectional  geometry  shown  as  first  column  in  Fig.  4(a))  used  in  the  experimental  studies. 


4.  Results 

Reconstructed  \la  and  jn's  images  obtained  from  the  noisy  simulated  data  with  imperfect  (7% 
error)  glandular  layer  priors  using  the  methods  described  in  Sec.  2  are  shown  in  Fig.  1.  Using 
hard  priors,  the  total  number  of  unknowns  is  reduced  to  6  parameters  (jla  and  /l'  for  each  of  the 
3  regions)  and  these  images  are  presented  in  the  last  column  of  Fig.  1(a).  The  images  from  two 
different  approaches  of  soft-priors  are  shown  in  the  middle  2  columns;  the  first  column  shows 
the  expected  results.  For  the  Helmholtz  case,  K  =  1/8  mm-1  (BPE)  was  used,  where  8  mm  is 
the  diameter  of  the  tumor.  Cross-sectional  plots  of  the  reconstructed  jda  and  /i'  distributions 
along  the  dotted  line  in  Fig.  1(a)  (first  column)  are  provided  in  Fig.  1(b).  Table  1(a)  and  (b) 
show  the  mean  and  standard  deviation  of  the  optical  property  estimations  in  each  region  of  the 
reconstructed  images.  Note  that  the  NIR  reconstruction  procedure  without  prior  information 
(not  shown)  did  not  generate  meaningful  images  in  this  complex  case. 

Figure  2(a)  shows  the  reconstructed  images  from  the  data  with  5%  noise  using  all  four  tech¬ 
niques  described  in  Sec.  2.  The  first  column  of  Fig.  1(a)  shows  the  actual  distribution  of  optical 
properties.  Figure  2(b)  shows  the  mean  and  standard  deviation  values  (as  error  bars)  of  recon¬ 
structed  images  using  different  techniques  for  different  regions  of  the  breast  geometry  with 
increasing  data  noise  level.  In  the  Helmholtz  case,  k  =  1/8  mm-1  (BPE)  was  used.  The  ac¬ 
tual  values  are  also  plotted  for  the  comparison.  It  can  be  clearly  seen  Laplacian  regularization 
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gives  lesser  standard  deviation  for  the  absorption  coefficient  ( fda )  reconstruction  compared  to 
the  Helmholtz  structure.  On  the  other  hand,  Helmholtz  structure  performs  better  than  Laplacian 
in  the  case  of  scattering  coefficient 

A  photograph  of  the  phantom  used  to  collect  data  at  785nm  with  16  fibers  in  a  single  plane 
(giving  240  measurements)  is  shown  in  Fig.  3.  Images  obtained  from  the  procedures  described 
in  Sec.  2  are  presented  in  Fig.  4(a)  along  with  cross-sectional  plots  of  the  reconstructed  ila 
and  jUy  distributions  in  Fig.  4(b).  Table  2(a)  and  (b)  show  the  mean  and  standard  deviation 
of  the  optical  property  estimations  in  each  region  of  the  reconstructed  images.  Here  also,  for 
Helmholtz  case,  the  BPE  (k  =  1/16  mm-1)  was  used.  Figure  5(a)  gives  the  results  for  different 
values  of  K  (given  on  the  top  of  the  each  column,  true  distribution  is  shown  as  first  column 
in  Fig.  4(a))  in  the  Helmholtz  case.  Corresponding  cross-sections  are  plotted  in  Fig.  5(b).  The 
mean  and  standard  deviation  values  for  each  of  this  case  are  also  given  in  Table  2(a)  and  (b). 

5.  Discussion 

The  reconstructed  results  (Fig.  1,  2,  4  and  5)  show  that  the  structural  priors  improve  the  re¬ 
constructed  image  quality  dramatically.  The  penalized  problem  formulation  (different  type  of 
regularizations)  generates  smoother  images  resulting  in  smaller  standard  deviations  from  the 
mean  values  (see  Table-2)  as  compared  to  the  generalized  problem  that  does  not  incorporate 
prior  information. 

The  hard-prior  case  produces  significant  optical  property  value  error  when  the  structural  a 
priori  information  is  imperfect  in  the  breast  geometry  (Fig.  1  and  Table- 1).  In  this  case,  a  7% 
variation  in  the  definition  of  the  glandular  layer  caused  false  estimates  of  the  optical  proper¬ 
ties.  On  the  other  hand,  soft-priors  (Laplacian  and  Helmholtz)  yield  good  estimates  for  each 
layer.  Hard-priors  over-estimate  the  tumor  absorption  coefficient  by  360%  and  under-estimate 
its  scattering  coefficient  by  88%  (Fig.  1).  Soft-priors  are  within  6%  of  the  expected  values  even 
with  the  error  in  the  structural  prior.  Note  that  in  this  particular  case,  hard-priors  failed  when 
7%  of  glandular  layer  made  as  a  fatty  layer,  yet  below  this  error  value  it  gave  reasonable  es¬ 
timates  of  optical  properties.  It  is  also  important  to  realize  that  this  is  only  a  test  case  of  the 
spatial-prior  error  progation.  There  is  not  enough  in-vivo  data  for  a  systematic  investiagion  of 
spatial-error  propagation  in  the  reconstruction  procedures,  especially  in  the  hard-priors  case. 

With  perfect  structural  priors,  increasing  the  data  noise  level  also  increases  the  quantification 
error  in  the  reconstructed  images  (Fig.  2(a)  and(b)).  Hard-Priors  clearly  fail  in  quantifying 
the  tumor  optical  properties  as  the  data  noise  level  increases.  Soft-Priors  does  better  in  the 
quantification  than  Hard-Priors.  It  is  also  evident  that  incorporation  of  structural  information 
is  key  for  accurate  quantification  of  the  optical  properties.  The  experimental  results  (Fig.  4) 
from  the  three-layer  gelatin  phantom  (Fig.  3)  also  show  that  incorporation  of  perfect  structural 
information  improves  the  quantification  and  quality  of  the  reconstructed  images.  Specifically, 
the  mean  and  standard  deviation  of  the  reconstructed  optical  property  values  in  each  region  are 
both  more  accurate  and  more  precise  where  the  priors  are  included.  The  mean  values  (Table-2) 
show  that  the  absorption  coefficient  for  the  Region-0  (fat)  is  under-estimated  by  77%  in  the  case 
of  the  Helmholtz  regularization  (BPE).  This  error  can  be  explained  by  the  fact  that  the  Euclidean 
distance  between  the  nodes  was  used  rather  than  the  distance  between  the  nodes  along  the 
boundary  of  a  particular  region.  Both  Laplacian  and  Hard-Priors  over  estimate  the  scattering 
coefficient  of  the  tumor  region  by  a  factor  of  2.  It  is  known  that  the  photon  path  length  is  affected 
by  the  scattering  coefficient.  By  constraining  the  problem  based  on  the  distance,  one  can  expect 
to  estimate  the  scattering  better.  In  this  experimental  case,  table-2  indicates  that  the  Helmholtz 
technique  produces  more  quantitative  accuracy  for  the  scattering  coefficient  estimation  of  the 
tumor  and  the  Laplacian  technique  is  best  for  the  absorption  coefficient  estimation  (as  well 
as  Fig.  2  and  4).  Usage  of  the  Helmholtz-type  form  for  the  diffusion  part  and  the  Laplacian- 
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Fig.  4.  (a)  Actual  fia  and  ii's  distributions  (axial  cross-section)  of  phantom  (Fig.  3)  case  are 
shown  in  the  first  column.  Optical  properties  for  the  region  labeled  ‘O’  are:  fia  =  0.0065 
mm-1  and  /i'  =  0.  65  mm-1.  Region  ‘1’  values  are:  fda  =  0.01  mm-1  and  fd's  =  1.0  mm-1. 
Region  ‘2’  (tumor)  values  are:  \ia  =  0.02  mm-1  and  fAfs  =  1.2  mm-1.  Reconstructed  jia  and 
ju'  distribution  from  different  techniques  (discussed  in  Sec.  2)  from  the  experimental  phan¬ 
tom  data.  Second  column  of  images  does  not  use  prior  information.  The  middle  rows  use 
soft  prior  structural  information  and  the  last  row  of  images  were  recovered  with  hard  priors. 
In  the  Helmholtz  case,  K  =  1/16  mm-1  (BPE)  was  used,  (b)  Cross-sectional  plots  along  the 
dotted  line  in  the  actual  image  (see  first  column  of  (a))  of  the  true  and  reconstructed  jj,a  and 
li’s  distributions. 
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Fig.  5.  (a)  Reconstructed  jj,a  and  /i'  images  from  the  experimental  phantom  data  using 
Helmholtz-type  regularization  matrix  for  different  values  of  K ,  which  are  given  at  the  top 
of  each  column,  (b)  Cross-sectional  plots  along  the  dotted  line  of  the  actual  images  in  Fig. 
4(a)  (first  column)  are  shown  with  the  data  from  reconstructed  fia  and  /i'  images  in  (a).  The 
best  prior  estimate  (BPE)  case  (k  =  1/16  mm-1)  is  also  presented  for  comparison. 
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type  form  for  the  absoprtion  part  (in  Eq.  8),  may  lead  to  a  better  estimation  of  both  optical 
properties,  which  will  be  considered  as  a  future  work.  Under  estimation  of  optical  properties 
in  the  regions  ‘O’  and  ‘V  can  be  attributed  to  the  random  fluctuations  in  heating/cooling  of 
the  gelatin  solution  especially  in  making  layered  phantom,  which  eventually  changes  optical 
properties.  Moreover,  ink  and  T1O2  particles  settling  to  the  bottom,  making  the  mixture  less 
homogenous  also  causes  the  change  optical  properties.  These  type  of  errors  are  more  common 
for  layered  gelatin  phantoms[37]. 


Table  2.  Mean  and  standard  deviation  of  the  reconstructed  (a)  jj,a  and  (b)  /i'  values  in 
different  regions  (labeled  in  first  column  of  Fig.  4(a))  recovered  from  the  experimental 
phantom  data.  The  corresponding  reconstructed  images  are  shown  in  Fig.  4(a)  and  5(a). 


Methods 

Region-0 

Region- 1 

Region-2 

Actual 

0.0065 

0.01 

0.02 

No  Priors 

0.0025±0.0010 

0.0045  ±0.0022 

0.0120±0.0090 

Laplacian 

0.0031  ±0.0002 

0.0051  ±0.0005 

0.0174±0.0029 

Helmholtz  (k  =  1/16) 

0.0015±0.0005 

0.0058±0.0009 

0.0241  ±0.0043 

Hard  Priors 

0.0032 

0.005 

0.0213 

Helmholtz  (k  =  1/5) 

0.0009±0.0006 

0.0061  ±0.0008 

0.0191±0.0031 

Helmholtz  (k  =  1/43) 

0.0027±0.0003 

0.0052±0.0007 

0.0234±0.0043 

Helmholtz  (k  =  1/86) 

0.0022±0.0005 

0.0061  ±0.0032 

0.0192±0.0044 

(a) 


Methods 

Region-0 

Region- 1 

Region-2 

Actual 

0.65 

1.0 

1.2 

No  Priors 

0.64±0.40 

0.66±0.27 

0.76±0.16 

Laplacian 

0.38±0.03 

0.63±0.07 

2.37±0.41 

Helmholtz  (k  =  1/16) 

0.46±0.02 

0.59±0.02 

1.08±0.12 

Hard  Priors 

0.37 

0.63 

2.74 

Helmholtz  (k  =  1/5) 

0.60±0.01 

0.57±0.01 

0.82±0.06 

Helmholtz  (k  =  1/43) 

0.39±0.03 

0.63±0.03 

1.19±0.13 

Helmholtz  (k  =  1/86) 

0.39±0.03 

0.63±0.04 

1.21±0.14 

(b) 


Theoretically  Helmholtz  and  Laplacian  structures  are  identical  when  K  =  0  (equivalently  / 
is  large).  Figures  5(a)  and  (b)  (as  well  as  Table-2)  show  that  when  k  =  1/86  mm-1  (/  =  86 
mm  is  the  diameter  of  the  domain),  Laplacian  and  Helmholtz  structures  give  reasonably  close 
reconstruction  values  of  optical  properties,  which  indicates  the  expected  trend  presented  in  this 
paper  for  the  two  methods.  It  also  indicates  that  the  BPE  case  (k  =  1/16  mm-1)  gives  the  best 
results,  as  the  priors  applied  are  close  to  the  true  structure  of  the  feature  (tumor).  These  results 
also  indicate  that  unreasonable  constraints  (like  K  =  1/5  mm-1 ,  Fig.  5(a)  first  column)  makes  the 
estimation  problem  amplify  the  noise  resulting  in  physiologically  invalid  (scattering  coefficient 
is  greater  for  fatty  layer  compared  to  fibroglandular  layer)  estimates  of  optical  properties.  In 
the  tumor  region,  as  k  decreases,  the  estimated  values  of  optical  parameters  become  closer  to 
expected  values  (Table-2  and  Fig.  5,  k  =  1/43  and  k  =  1/86  cases).  This  is  due  to  the  correlation 
length  becoming  larger,  making  the  covariance  in  the  neighboring  nodes  larger. 

In  this  study,  it  is  shown  that  imperfect  priors  (commonly  caused  by  improper  image  seg¬ 
mentation  and  image  artifacts  in  MRI  or  X-ray)  can  lead  to  error-prone  results  in  the  hard-prior 
case  whereas  soft-priors  are  more  immune  to  uncertainty  in  the  prior  data.  It  is  also  shown  that 
the  techniques  used  to  incorporate  the  soft  structural  prior  information  influences  the  image 
outcome,  which  may  lead  to  improvements  in  image  accuracy  if  properly  implemented.  Srini- 
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vasan  et.  al[38]  have  found  that  5%  error  in  optical  properties  introduces  approximately  45% 
error  in  the  chromophore  concentration  when  a  limited  number  of  wavelengths  are  used  for 
imaging.  The  correct  “soft”  inclusion  of  a  priori  information  therefore  can  be  expected  to  lead 
to  a  more  accurate  quantification  of  chromophore  concentrations  as  well.  It  should  be  noted  that 
over- weighting  of  the  penalty  term  in  the  problem  formulation  may  make  the  solution  ignore 
the  data-model  misfit  and  emphasize  smooth  feature  extraction.  The  techniques  developed  in 
this  work  were  applied  for  two-dimensional  test  objects,  and  can  be  easily  extended  to  three- 
dimensional  case.  A  more  extensive  study  of  this  is  left  for  future  investigations. 

6.  Conclusions 

This  work  has  investigated  several  ways  of  incorporating  structural  information  into  an  itera¬ 
tive  image  reconstruction.  The  results  have  been  supported  by  gelatin  phantom  experiments  that 
represent  multi-layered  structure  which  is  commonly  found  in  breast  tissue  with  adipose  (fatty) 
tissue  on  the  exterior  and  fibroglandular  tissue  nearer  to  the  interior.  Soft-priors  allow  the  tissue 
optical  properties  to  vary  within  predefined  regions,  unlike  hard-priors  which  constrain  these 
zones  to  be  homogeneous.  Hard-priors  were  found  to  perform  poorly  when  the  prior  informa¬ 
tion  contained  area  errors  as  small  as  7%  which  can  easily  be  produced  by  most  segmentation 
algorithms.  True  boundary  extraction  from  MRI  images  introduces  unavoidable  segmentation 
and  discretization  errors  that  are  better  tolerated  when  the  structural  information  is  encoded 
through  the  soft-prior  approach  involving  a  penalty  matrix. 

The  results  reported  here  indicate  that  the  optical  properties  of  different  tissue  types  can 
be  quantified  more  accurately  when  their  estimation  is  properly  guided  by  ’’soft”  structural 
a  priori  information.  The  problem  formulation  and  results  presented  in  this  work  indicate 
that  data  from  other  imaging  modalities  such  as  ultrasound  or  x-ray  tomography,  could 
also  be  used  as  the  source  of  the  structural  prior.  In  the  experimental  cases  investigated, 
the  Helmholtz  structure  using  best  priors  gave  a  better  estimation  of  scattering  coefficient 
of  the  target  (tumor).  However,  the  Laplacian  type  of  regularization  leads  to  more  superior 
absorption  coefficient  estimate.  So  it  is  reasonable  to  conclude  that  Laplacian  structure  gives 
the  best  estimates  of  total  hemoglobin  concentration  ( Hbj ),  hemoglobin  oxygen  saturation 
(St 02%)  and  water  fraction  ( H2O )  (which  are  the  main  absorbers).  Helmholtz  structure  gives 
the  best  estimates  of  the  scattering  power  and  scattering  amplitude  (scattering  parameters). 
The  framework  presented  here  can  also  be  extended  to  other  regularization  terms  such  as  total 
variation  minimization  or  spectral  prior  constraints,  which  may  be  studied  in  future  work. 
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ABSTRACT 


Three-dimensional  (3D)  diffuse  optical  tomography  (DOT)  is  known  to  be 
a  non-linear,  ill-posed  and  sometimes  under-determined  problem,  where 
regularization  is  added  to  the  minimization  to  allow  convergence  to  a  unique 
solution.  In  this  work,  a  generalized  least-squares  (GLS)  minimization  method 
was  implemented,  which  employs  weight  matrices  for  both  data-model  misfit 
and  optical  properties  to  include  their  variances  and  covariances,  using  a 
computationally  efficient  scheme.  This  allows  inversion  of  a  matrix  that  is 
of  dimension  dictated  by  the  number  of  measurements,  instead  of  by  the 
number  of  imaging  parameters.  This  increases  the  computation  speed  up  to 
four  times  per  iteration  in  most  of  under-determined  3D  imaging  problems. 

An  analytic  derivation,  using  the  Sherman-Morrison- Woodbury  identity,  is 
shown  for  this  efficient  alternative  form  and  it  is  proven  to  be  equivalent, 
not  only  analytically,  but  also  numerically.  Equivalent  alternative  forms  for 
other  minimization  methods,  like  Levenberg-Marquardt  (LM)  and  Tikhonov, 
are  also  derived.  3D  reconstruction  results  indicate  that  the  poor  recovery  of 
quantitatively  accurate  values  in  3D  optical  images  can  also  be  a  characteristic 
of  the  reconstruction  algorithm,  along  with  the  target  size.  Interestingly,  usage 
of  GLS  reconstruction  methods  reduces  error  in  the  periphery  of  the  image, 
as  expected,  and  improves  by  20%  the  ability  to  quantify  local  interior  regions 
in  terms  of  the  recovered  optical  contrast,  as  compared  to  LM  methods. 
Characterization  of  detector  (PMTs)  noise  have  enabled  the  use  of  the  GLS 
method  for  reconstructing  experimental  data  and  showed  a  promise  for  better 
quantification  of  target  in  3D  optical  imaging.  Use  of  these  new  alternative 
forms  becomes  effective  when  the  ratio  of  the  number  of  imaging  property 
parameters  exceeds  the  number  of  measurements  by  a  factor  greater  than  2. 

Keywords:  near  infrared,  diffuse  optical  tomography,  three-dimensional  imaging,  image  re¬ 
construction,  inverse  problems,  least  squares  minimization. 
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1.  Introduction 


Diffuse  optical  tomography  (DOT)  uses  near  infrared  wavelengths  (600-1000  nm)  to  obtain 
optical  absorption  and  scattering  images  for  characterizing  functional  properties  of  the 
tissue  under  investigation.1-4  The  most  important  step  in  forming  these  images  is  solving 
the  inverse  problem,  i.e.  estimating  the  optical  properties  by  matching  the  experimental 
data  with  modeled  results  in  the  least-squares  sense.3,5,6  This  problem  is  typically  ill-posed 
and  ill-determined  depending  on  the  noise  in  the  data,  the  number  of  measurements  and 
the  dimensions  of  the  parameter  space.6  Even  though  light  travels  in  three-dimensions  (3D), 
most  of  the  numerical  models  reported  in  the  literature  have  been  two-dimensional  (2D) 
because  of  computational  considerations.  Moreover,  the  3D  DOT  imaging  problem  is  more 
under-determined  relative  to  the  2D  case  and  has  been  found  to  generate  poor  quantitative 
estimates  of  the  optical  properties  when  compared  to  2D  results.7-14  Several  methods 
have  appeared  in  the  literature  describing  efficient  3D  computations,11,14, 16  but  no  unified 
approach  to  the  problem  has  been  discussed.  Recently,  a  generalized  least-squares  (GLS) 
minimization  scheme  was  presented  for  2D  DOT  image  reconstruction.17  This  paper  reports 
a  computationally  efficient  approach  for  implementing  GLS  minimization  in  3D  which  shows 
an  improvement  in  the  quantification  of  optical  properties  relative  to  earlier  studies.  Even 
though  the  focus  is  on  GLS  implementation,  equivalent  forms  for  other  methods  are  also 
presented  in  light  of  the  GLS  framework. 

The  inverse  problem  in  DOT  is  solved  by  minimizing  the  objective  function  (0)  over 
the  range  of  optical  properties  (/i).  Methods  based  on  gradient  optimization,  which  do  not 
require  an  explicit  inversion  of  the  Hessian  matrix  1  (in  general  some  form  of  JTJ,  J  being 
the  Jacobian)  are  known  to  be  computationally  efficient.18,19  But  these  methods  require 
an  optimization  scheme,  which  can  be  thought  of  as  an  inner  iteration,  for  choosing  the 
step-size,  and  are  not  as  straight  forward  as  a  direct  inversion  of  the  matrix.  Alternatively, 
the  full-Newton  methods  require  calculation  of  the  Jacobian  (J),  the  forward  data,  and 
inversion  of  the  dense  Hessian  matrix  at  each  iteration.  Because  full-Newton  methods  are 
relatively  easy  to  implement,  they  are  widely  used  for  DOT  image  reconstruction  even 
though  they  require  large  matrix  inversions  at  every  iteration.  Thus,  while  the  full-Newton 
method  is  ideal  for  small  problems,  it  rapidly  becomes  intractable  for  larger  domains,  such 
as  those  encountered  in  3D  imaging  problems.  This  manuscript  presents  a  formal  approach, 
using  the  Sherman-Morrison- Woodbury  identity,  to  construct  a  more  efficient  alternative 
forms  of  update  equations  of  GLS  and  Levenberg-Marquardt  (LM)  minimization  schemes, 
which  generates  a  Hessian  matrix  to  invert.  The  dimension  of  this  Hessian  is  dictated  by 
the  number  of  measurements,  rather  than  the  number  of  parameter  estimates  which  can 

1Here,  Hessian  approximates  the  second  derivative,  this  can  be  some  form  equivalent  to  JrJ  or  JJr 
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be  considerably  lower  for  highly  under-determined  problems,  and  therefore,  much  more 
efficient  computationally.  This  equivalent  form  is  also  shown  for  other  common  minimization 
method,  namely  Tikhonov. 

Later  part  of  this  paper  describes  a  way  to  characterize  the  systematic  noise  using  a  simple 
analytical  formula,  when  photomultiplier  tube  (PMT)  used  as  a  detector.  Characterizing 
noise  behavior  of  the  experimental  data  lead  to  use  of  the  GLS  technique  in  the  experimental 
data  case  and  it  is  also  shown  that  usage  of  noise  characterstics  will  lead  to  better  quality 
and  quantification  of  target  in  a  experimental  test  case. 

2.  Diffuse  Optical  Tomography:  Forward  Problem 

Near- infrared  (NIR)  light  propagation  in  a  biological  tissue  like  breast  can  be  modeled  using 
the  Diffusion  Equation  (DE)6,20  which  in  the  frequency  domain,  becomes 

— V.-D(r)V<h(r,  u)  +  (/ra(r)  +  iu/c)$(r,u)  =  Qa(r,uj )  (1) 

where  the  optical  diffusion  and  absorption  coefficients  are  given  by  D{r )  and  //„  (r) ,  respec¬ 
tively.  The  light  source,  represented  by  Q0{r,  u),  is  modeled  as  isotropic.  <f>(r,  u)  is  the  photon 
fluence  rate  at  a  given  position  r.  The  light  modulation  frequency  is  denoted  by  u,  where 
u  =  27rf,  (here  f  =  100  MHz).  The  velocity  of  light  in  tissue  is  represented  by  c,  which  is 
assumed  to  be  constant.  Note  that 


3  \j*a(r)  +  fJ.'a(r)] 

where  p/s(r)  is  the  reduced  scattering  coefficient  which  is  defined  as  p/s  =  ys(l  —  g).  /q  is 
the  scattering  coefficient  and  g  is  the  anisotropy  factor.  The  finite  element  method  (FEM)  is 
used  to  solve  equation  (1)  to  generate  modeled  data  (G(/r))  for  a  given  distribution  of  optical 
properties  (/i), 11,13,21  where  fi=  [fj,'s{r)]  na{r)\.  A  Type-III  boundary  condition  is  employed 
to  account  for  the  refractive- index  mismatch  at  the  boundary.22  Under  the  Rytov  approxi¬ 
mation,  the  data  (y)  becomes  the  natural  logarithm  of  the  amplitude  (A)  and  phase  ( 9 )  of 
the  frequency-domain  signal;  y  =  [lnA;6\. 

3.  Diffuse  Optical  Tomography:  Inverse  Problem 

3. A.  Levenberg-Marquardt  (LM)  Minimization 

The  most-common  approach  for  solving  the  inverse  problem  in  DOT  is  Levenberg-Marquardt 
(LM)  minimization.1,6,11,13,17,20,23,24  A  detailed  discussion  of  this  method  is  available  in 
Ref.17  and  it  is  only  briefly  reviewed  here. 
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The  objective  function25, 26  for  this  approach  is  defined  as 


n  =  \\y  -  G{y)f 


(3) 


where  y  is  the  experimental  data  and  G(/i)  is  the  modeled  response.  Minimization  of  this 
objective  function  with  respect  to  p  is  achieved  by  setting  the  first-order  derivative  equal  to 
zero 


^  =  3tS  =  0 

Ofl 


(4) 


where  8  is  the  data-model  misfit,  8  =  y  —  Gi/i),  and  J  represents  the  Jacobian  (J  =  - ) . 
Due  to  the  ill-conditioned  nature  of  the  problem,  the  update  equation  for  the  optical  prop¬ 
erties  at  iteration  ‘i’  is  regularized  to  be 


A  &  =  [JTJ  +  al]  1  JTSi_1 


(5) 


or  equivalently  (See  Appendix- A. 2) 

Apj  =  JT  [JJT  +  al] -1  Si-i  (6) 

where  A//,.,  represents  the  update  of  the  optical  property  parameters  at  the  ith  step,  a  is  the 
regularization  parameter,  which  monotonically  decreases  with  increasing  iteration  (always 
>  0).26  In  this  approach  (Eq.  5),  the  Jacobian  is  normalized  by  its  optical  properties. 
Moreover,  a  is  chosen  empirically  (it  typically  starts  at  10  and  reduced  by  a  factor  of  10°'25 
at  every  following  iteration  after  being  multiplied  by  the  maximum  of  the  diagonal  values  of 
JJT.17,27).  The  iterative  procedure  is  stopped  when  the  L2  norm  of  the  data-model  misfit 
(J)  does  not  improve  (in  our  experience,  by  more  than  1%  because  beyond  these  values  the 
LM  procedure  can  become  unstable17). 


Even  though  LM  minimization  or  its  modified  versions  have  been  used  successfully  for 
DOT  image  reconstruction,1,6,11,13,17,20,23,24,27  the  final  image  depends  on  the  choice  of  a 
due  to  the  ill-conditioned  nature  of  the  problem.  Moreover,  the  approach  ignores  the  noise 
characteristics  of  the  data  and  optical  properties.  A  more  systematic  and  generalized  method 
for  image  reconstruction  can  be  based  on  GLS  minimization.  The  GLS  scheme  is  discussed 
extensively  in  Ref.17  and  is  only  briefly  reviewed  here. 

3.B.  Generalized  Least  Squares  (GLS)  Minimization 
In  GLS,  the  objective  function  is  given  by17,28 

Ll  =  (y-  G(n))TWs{y  -  G(fj))  +  (p  -  /a 0)TW#1_M (p  -  /i0)  (7) 

where  W,j  and  W/t_/i0  are  the  weight  matrices  for  the  data-model  misfit  (J)  and  optical 
properties  (/.i-//o),  respectively.  Note  that  W $  =  (C^)-1,  where  C  represents  the  covariance 
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matrix,  and  similarly  (see  Appendix: A-4  of  Ref.28).  These  weight  ma¬ 

trices  are  symmetric  and  positive  definite  (because  they  are  inverses  of  covariance  matrices). 
No  regularization  parameter  is  involved  because  the  weight  matrices  include  the  noise  char¬ 
acteristics  of  the  experimental  data  and  optical  properties.17  Similarly  to  the  LM  approach, 
minimization  of  fl  (Eq.  7)  is  accomplished  by  setting  the  first  derivative  of  0  with  respect 
to  n  equal  to  zero: 

BO 

—  =  JT W s5  -  W^0  {n  -  no)  =  0.  (8) 

Linearizing  the  problem  leads  to  the  iterative  update  equation  (for  ith  iteration)17 

A m  =  [JTW5J  +  (JTW5h,_!  -  -  li o))  (9) 

Explicit  definitions  of  the  weight  matrices  (Wj  and  are  given  in  Ref.17  Although 

any  number  of  forms  for  'Wli-I,0  can  exist,  only  one  is  considered  here,  specifically,  where 
the  covariance  matrix  is  defined  as17 

[CM-Mo]q  =  K-m)2  (i  +  y)  e“ ^  (10) 

with  l  being  the  correlation  length  (here  l  =  15  mm)  and  r,tj  being  the  distance  between  the 
FEM  nodes  i  and  j.  (cr/f_//0)2  is  the  expected  variance  of  //  —  /if,-  Strategies  for  calculating 
the  expected  variance  are  given  in  Ref.17  In  this  work,  the  expected  variance  is  determined 
from  the  prior  knowledge  that  the  expected  contrast  between  tumor  and  normal  tissue  is 
about  50-400%.  To  demonstrate  the  robustness  of  the  GLS  reconstruction  procedure,  for  the 
results  discussed  here,  the  variance  was  chosen  to  be  (4*  /i )2.  Both  weight  matrices,  W()-  and 
WM_M0,  are  computed  before  the  reconstruction  procedure  begins,  whereas  the  Jacobian  (J), 
and  modeled  data,  are  calculated  at  each  iteration.  The  iterative  procedure  is  stopped 

when  the  L2  norm  of  the  data-model  misfit  (J)  does  not  improve  by  more  than  0.001%. 
Beyond  these  values,  the  round-off  error  dominates. 

3.B.I.  GLS  Implementation 

The  parameters  recovered  in  the  case  of  this  GLS  scheme  are  (//(,;  fj,a) ,  which  is  different 
from  some  previous  approaches  that  estimate  (D;//a).  The  later  case  has  a  mismatch 
because  the  units  of  D  are  mm  whereas  those  of  //n  are  mm-1.  In  its  implementation 
typically  the  whole  equation  is  normalized  by  the  optical  properties  (outlined  in  Ref.17) 
which  becomes  computationally  intensive  especially  for  GLS  in  3D  because  the  update 
equation  must  be  left  and  right  multiplied  by  the  optical  properties  at  every  iteration. 
Here,  the  GLS  problem  was  reformulated  in  terms  of  (/r's;  //a ) ,  so  that  both  parameters  have 
same  units  (mm-1).  While  this  is  a  relatively  minor  alteration  in  the  form  of  the  algorithm, 
it  has  important  implications  for  the  computational  time  required  for  matrix  preconditioning. 
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A  simple  transformation  converts  the  diffusion  part  of  the  Jacobian  (Ajfp-)  to  its  scattering 
component  ( 9g^ ) : 

dG{n)  dG(n)  dD 


Using  Eq.  2 


dl-i's 

dD  1 


dD  dn's 


-1 


=  - 3D 2 


dfi's  3  \[(fia  +  /i')]2 
and  substituting  Eq.  12  in  Eq.  11,  leads  to 

dGjjj)  _  dGjjj ’,)  2 

dfl's  dD  (  j 

After  this  transformation  (Eq.  13)  the  Jacobian  (J)  has  the  form 


dlnAi  dlnAi 

dlnAi 

dlnAi 

dlnAi 

dlnAi 

3/4 1  dn's2 

d^'sNN 

dual 

d^a2 

dflaNN 

d9i  dji 

<90i 

<90q 

<90! 

dQi 

■Vsi  &n's2 

dv'sNN 

d[lal 

d[la2 

d^aNN 

dlnA<2  dlnA2 

dlnA2 

dlnA2 

dlnA2 

dlnA2 

■V3i  dv',2 

d^'sNN 

d/iai 

dna2 

d^aNN 

<902  <902 

<902 

d02 

d02 

<902 

VS1  dv's2 

d^'sNN 

dflal 

d^a2 

dflaNN 

din  An  m  din  An  m 

din  An  m 

din  An  m 

din  An  m 

din  An  m 

■VS1  dv's2 

d^'sNN 

d^al 

dna2 

dflaNN 

99nm  99nm 

<VS  1  0/ 42 

dON  m 

d^'sNN 

90nm 

dON  m 

dON  m 

dual 

d{la2 

d^aNN 

(11) 

(12) 

(13) 


(14) 


where  NM  and  NN  represents  the  number  of  measurements  and  number  of  property  param¬ 
eters  associated  with  the  FEM  mesh,  respectively.  The  implementation  of  the  GLS  update 
equation  (Eq.  9)  requires  assembly  of  the  weight  matrix  (W/(_/icl )  for  simultaneous  recon¬ 
struction  of  fj/s  and  fj,a  and  is  accomplished  by  writing  Eq.  9  in  block  matrix  form 


H< 

H v's^a 

4_ 

w  ,  , 

0 

\ 

A/Lig  ^ 

(jt: 

)^W55i_1 

w ^'s-^'so  Mso) 

Hu  2 

Ma  J 

\ 

0 

^a-MaO 

J 

_  uT: 

W  W85i-1 

Wpa-MaO  (^ai- 1  —  A^ao) 

(15) 

where  H  represents  the  Hessian  matrix  ( JTWf,J) .  Here,  the  cross-terms  in  the  weight 
matrix  (W^J  are  zero  because  /r'  and  //,,  are  independent  parameters  in  the  estimation 
procedure.  Note  that  the  dimensions  of  the  matrices  in  Eq.  9  are:  J:  2NM  x  2NN,  WM_/(0 : 
2NN  x  2NN,  Wj:  2NM  x  2NM,  8  =  2NM  x  1,  and  A/i  =  2NN  x  1.  Most  3D-DOT  problems 
are  ill- determined,  i.e.  NM  <C  NN. 


Computing  an  update  of  the  optical  properties  (A/q,  from  Eq.  9  )  requires  an  inversion 
(or  its  equivalent)  of  a  large  matrix  with  dimensions  2NN  x  2NN.  Inverting  a  matrix  of 
dimension  NxN,  typically  requires  an  order  of  N3  operations  and  N2  memory.28  Hence,  any 
gain  in  reducing  the  dimensionality  of  the  matrix  to  be  inverted,  will  reduce  the  computation 
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time  cubically  and  the  memory  requirement  quadratically.  An  alternative  from  of  Eq.  9,  which 
requires  few  operations  is 


A m=  i-c/1_/J0jt(jc ^0jt  +  c5)  1 J  (16) 


Full  derivation  of  this  alternative  form  is  given  in  the  appendix  (A.l)  along  with  the 
equivalent  expressions  for  other  minimization  methods.  Eq.  16  requires  an  inversion  of  a 
matrix  with  dimensions  2NM  x  2NM  (same  is  true  for  Eq.  6). 


Note  that  the  covariance/weight  matrices  are  calculated  before  the  start  of  the  iterative 
procedure  and  are  used  throughout  the  iteration.  For  nodes  where  the  sensitivity  (column 
sum  of  Jacobian)  fell  below  1%  of  the  maximum  sensitivity,  the  expected  variance  of  the 
optical  properties  was  chosen  to  be  1%  of  background  fi  (in  Eq.  10). 

4.  Simulation  Studies:  Three-dimensional  Test  Problem 

For  all  numerical  experiments  discussed  here,  the  imaging  domain  was  chosen  to  be 
cylindrical  (as  shown  in  Fig.  1)  with  a  diameter  of  84  mm  and  height  of  109  mm.  The 
background  optical  properties  were  fj,a  =  0.01  mm-1  and  /j/  =1.0  mm-1.  Two  meshes  were 
used:  (1)  a  cylinder  consisting  of  21,440  nodes  corresponding  to  110,483  linear  tetrahedral 
elements  for  the  forward  model  and  (2)  a  cylinder  having  9,211  nodes  corresponding  to 
45,980  linear  tetrahedral  elements  for  the  reconstruction.  The  data-collection  geometry 
consisted  of  48  fibers  that  were  arranged  in  a  circular,  equally-spaced  fashion  in  three  layers 
spaced  10  mm  apart  (Fig.  1),  with  16  fibers  per  plane.  One  fiber  was  used  at  a  time  as  the 
source  while  the  fibers  in  the  same  “source  fiber  plane”  were  used  as  detectors  to  generate 
720  (3x16x15)  measurement  locations  or  a  total  of  1440  values  (720  InA  data  points  and 
720  9  data  points).  The  sources  were  modeled  as  a  Gaussian  profile  with  a  full- width  half 
maximum  of  3  mm  to  represent  the  distribution  used  in  an  experimental  setup.29  The  source 
was  also  placed  one  mean  transport  scattering  distance  inside  the  boundary. 


Both  spherical  and  cylindrical  objects  were  considered  as  targets.  The  cylindrical  target 
had  a  contrast  of  2:1  with  the  background  in  both  fj,a  and  //(  and  a  diameter  of  15  mm.  It 
extended  in  Z-direction  throughout  the  domain  (height  of  109  mm)  and  was  placed  at  the 
center  (at  (0,0),  first  row  of  Fig.  5)  and  near  the  boundary  (at  (30,0)).  Two-dimensional 
cross-sections  of  both  reconstructed  and  actual  3D  volumes  are  displayed  in  increments  of 
5  mm  spanning  from  z  =  -25  mm  to  z  =  25  mm  (from  left-hand  side  to  right-hand  side) 
are  shown  in  Figs.  3,  5,  and  6.  Cross-section  below  z=  -25  mm  and  beyond  z  =  25  mm 
did  not  show  any  deviation  from  the  starting  values  of  iterations  as  the  sensitivity  in  this 
region  is  almost  negligible  compared  to  the  rest  of  the  domain,  so  these  cross-sections  are 


omitted  for  display  purposes.  Measurements  with  a  noise  level  of  1%  were  assumed  as  the 
experimental  data  (y)  in  most  of  the  cases  discussed  here.  The  noise  variance  was  also 
used  with  GLS  reconstruction  algorithm.17  The  background  optical  properties  were  selected 
as  starting  values  for  the  iterative  image  reconstruction  procedures  discussed  in  Sec. 3.  All 
computations  were  carried  out  on  a  Linux  work  station  with  an  AMD  Dual  Core  Opteron 
280  processor  (2.2  GHz)  with  8GB  of  RAM. 

5.  Phantom  Studies 

5. A.  Data  variance  estimation 

Use  of  weight  matrices  in  the  GLS  scheme  (Wj  in  Eq.  9)  requires  an  estimation  of  data 
variance,  which  requires  experimental  characterization  of  the  expected  values,  prior  to 
patient/phantom  imaging.  This  was  achieved  by  tracking  the  detected  voltage  measured 
at  the  photo  multiplier  tube  (PMT).  PMTs  are  used  as  a  detectors  in  the  experimental 
system  at  Dartmouth,  details  of  the  experimental  system  are  given  in  Ref.29  Note  that  this 
characterization  includes  only  systematic  errors  associated  with  low  signal  levels,  but  errors 
due  to  poor  fiber-tissue  coupling  are  not  accounted  for  in  this  model. 

Starting  from  the  assumption  that  the  detected  signal  using  a  PMT  in  diffuse  optical 
imaging  is  shot  noise  limited  leads  to 


a  =  VN  (17) 

where  a  is  the  standard  deviation  in  the  detected  signal  and  N  is  the  number  of  photons 
reaching  the  PMT.  The  voltage  (V)  measured  at  the  PMT  is  directly  proportional  to  N, 
which  also  implies  that  amplitude  ( A )  of  the  detected  frequency  domain  signal  (y)  is  pro¬ 
portional  to  this  voltage.  This  is  written  as 

A  ocV  =>  A  =  kV 

(18 

cr(A)  =  ka(V) 

here,  k  acts  as  a  proportionality  constant.  In  the  reconstruction  procedure,  the  Rytov  ap¬ 
proximation  is  used,  leading  to  data  being  represented  as  InA  rather  than  A.  If  f(x)  is  a 
function  of  x  and  is  continuous  and  differentiable,  then 

°(f(x))  =  m^<r(x)  (19) 

similar  to  the  previous  equation  (Eq.  19),  writing  the  standard  deviation  of  InA  leads  to 

a  (InA )  =  —a(A)  (20) 
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now  using  Eq.  18  leads  to 


a(lnA)  =  (21) 

From  the  above  equation  (Eq.  21),  it  can  be  concluded  that  the  variance  in  data  (a2)  can 
be  known  by  measuring  the  deviation  in  the  PMT  voltage  (V). 


To  measure  the  deviation  in  the  measured  signal,  a  series  of  light  signal  measurements 
were  taken  through  a  homogeneous  intralipid  phantom  experiments  were  conducted  with 
increasing  levels  of  blood  (HbT)  concentration,  varying  from  7.3  /iM  to  36  /iM,  leading 
to  a  decrease  in  the  measured  PMT  voltage.  To  achieve  this,  the  gain  of  the  PMT  was 
kept  at  0.9.  A  concise  discussion  of  the  PMT  gain  setting  in  the  system  characterization  is 
given  in  Ref.29  A  single  source  and  the  farthest  detector  was  used  for  these  transmission 
measurements.  For  every  concentration,  200  data  points  were  collected  to  estimate  the 
deviation  in  the  measured  voltage  using  the  same  gain  settings.  The  approach  for  the 
characterization  is  similar  to  the  one  described  in  Ref.,29  except  the  raw  detected  voltage 
was  used  here  for  estimation  of  the  error  (or  deviation  a).  Note  also  that  two  sets  of 
diameters,  56  mm  and  84  mm,  were  used  to  get  the  voltage  in  the  range  of  0-1  volts.  This 
was  repeated  for  all  the  wavelengths  to  ensure  uniformity  of  performance  in  the  signal,  and 
to  ensure  that  the  obeserved  trend  was  independent  of  wavelength  and  gain  setting. 


Figure  9  gives  a  plot  of  error  (cr(V)/V)  as  a  function  of  measured  PMT  voltage  for  785 
nm  wavelength.  A  similar  trend  was  observed  for  other  wavelengths.  This  plot  also  gives  a 
deviation  in  phase  (cr(0))  in  degrees  for  the  same  voltage.  Each  of  these  points  represents 
a  sample  size  of  200.  The  lowest  measured  voltage  of  PMT  was  0.001  volts.  The  measured 
deviations  were  1%  for  InA  and  0.5°  for  6  for  PMT  voltages  above  0.005  volts.  These  values 
are  similar  to  the  ones  reported  in  the  literature  (Mcbride  et.  al29  reported  0.32%  for  the 
PMT  voltage  and  0.48°  in  phase).  A  solid  line  in  Fig.  9  shows  these  average  deviation  using 
1/V2  fitting  (following  the  shot-noise  model).  From  this  plot,  the  weight  matrix  for  the 
data- model  misfit  (Wj  in  Eq.  9)  can  be  written  as 


[W,],,  = 


/ 

0 

if  i 

%  j 

[(5.7/V2)+0.6]2  f°V 

InA  with 

V  >  0.001 

if  i 

=  j 

<  [(1.3/V2)+0.4]2  f°V 

9  with  V 

>  0.001 

if  i 

=  j 

1 

for  InA 

with  V  < 

0.001 

if  i 

=  j 

1 

l  W2 

for  6 

with  V  < 

0.001 

if  i 

=  j 

(22) 


This  implies  that  if  the  PMT  voltage  is  below  0.001  volts,  the  signal  is  considered  to  be  in 
the  noise  floor.  For  the  signals  above  0.001  volts,  the  variance  can  be  estimated  using  1/V2 
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fitted  curve  (Eq.  22).  This  characterization  enables  the  usage  of  the  GLS  scheme  in  the  case 
of  experimental  data,  with  a  requirement  that  the  voltage  measured  at  the  PMT  is  available. 

5. B.  Gelatin  Phantom 

A  multi-layer  cylindrical  gelatin  phantom  of  diameter  86  mm,  height  25  mm  was  fabricated 
using  different  mixtures  of  India  ink  for  absoprtion  and  Titanium  oxide  (TiC^)  for  scattering. 
These  different  layers  of  gelatin  were  fabricated  by  successively  hardening  heated  gelatin 
solutions  (typically  80%  deioinized  water  and  20%  gelatin  (G2625,  Sigma  Inc))  along  with 
different  amounts  of  ink  and  Ti02  (Sigma  Inc).  A  cylindrical  hole  extending  in  Z-direction 
(diameter  of  16  mm  and  height  of  24  mm)  filled  with  intralipid  mixed  with  india  ink  acted 
as  a  target  having  the  optical  properties  jia  =  0.02  mm-1  and  /i's  =  1.2  mm-1.  The  outer 
layer  with  optical  properties  pa  =  0.0065  mm-1  and  n's  =  0.65  mm-1  had  a  thinkness  of  10 
mm  mimicing  the  typical  fatty  layer  of  the  breast.46  The  middle  layer  with  76  mm  diameter, 
mimicing  fibroglandular  layer,  had  optical  properties  /ia  =  0.01  mm-1  and  /%  =  1.0  mm-1. 
Validation  of  individual  layers  optical  properties  was  performed  by  the  data  collected  on 
large  cylindrical  samples  of  each  layer  using  785  nm  wavelength  laser  diode  as  the  source. 
Two-dimensional  cross-sections  of  this  gelatin  phantom  optical  properties  are  displayed  in 
increments  of  2.5  mm  spanning  from  z  =  -12.5  mm  to  z  =  12.5  mm  (from  left-hand  side 
to  right-hand  side)  in  top  rows  of  Fig.  10  (a)  and  (b).  In  this  phantom  case,  data  was 
collected  using  only  one  layer  of  fibers  (at  z  =  0  mm)  leading  to  240  InA  data  points  and 
240  9  data  points.  A  cylindrical  mesh  consisting  8990  nodes  corresponding  to  44803  linear 
terahedral  elements  was  used  and  the  experimental  data  was  also  calibrated  using  a  reference 
homogenous  phantom  data.  The  outer  layer  (mimicing  fatty  layer)  optical  properties  were 
used  as  initial  guess  for  reconstruction  procedures.  A  second  mesh  with  the  same  geometry 
contianing  3718  nodes  (16627  linear  tetrahedral  elements)  was  used  as  a  reconstruction  mesh. 

6.  Results 

The  number  of  operations  required  to  produce  an  optical  property  update  (A/p)  for  both  the 
original  and  alternative  GLS  update  equations  (Eqs.  9  and  16,  respectively)  was  compared 
as  a  function  of  the  ratio  of  the  number  of  estimation  parameters  (N N)  to  the  number  of 
measurements  (NM).  A  similar  comparison  for  LM  update  equations  (Eqs.  5  and  6)  was  also 
performed.  The  results  are  plotted  in  Fig.  2(a)  and  (b)  respectively.  The  expressions  used 
for  calculating  the  number  of  operations  are  given  in  Appendix-A.4.  Memory  required  for 
implementing  these  inversions  are  presented  in  Fig.  2(c).  To  relate  the  analysis  to  existing 
experiments,  NM  was  chosen  to  be  720.  The  number  of  nodes  (equivalently  estimation 
parameters)  was  varied  from  2  to  80,000.  The  number  of  computations  increases  for  both 
LM  and  GLS  cases  as  NN  increases,  but  the  alternative  form  for  GLS  (GLS-AF)  has  a 
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lower  computational  cost  when  NN/NM  is  greater  than  6.  In  the  case  of  alternative  form 
LM  (LM-AF),  this  is  when  NN/NM  is  greater  than  2.  In  terms  of  memory,  as  soon  as 
NN  >  NM,  both  GLS-AF  and  LM-AF  requirement  is  less  than  GLS  and  LM  counterparts. 

In  order  to  assess  algorithm  performance,  a  series  of  test  reconstructions  were  evaluated. 
Data  with  1%  noise  from  the  cylinder  containing  a  15  mm  diameter  spherical  target  (first 
rows  of  Fig.  3(a)  and  (b))  was  reconstructed  for  the  optical  properties  using  LM  and 
GLS  techniques  (middle  and  last  rows  and  Fig.  3(a)  and  (b)).  In  the  case  of  the  GLS 
scheme,  the  reformulated  update  equation  (alternative  form  -  Eq.  16)  was  used  (last  rows 
of  Fig.  3(a)  and  (b)).  However,  it  was  also  important  to  confirm  that  the  two  forms  of 
the  update  equation  (Eqs.  9  and  16)  produced  numerically  equivalent  solutions.  Figure 
4  shows  a  comparison  of  results  generated  with  the  original  GLS  update  equation  (Eq. 
9)  and  its  alternative  form  (Eq.  16)  in  terms  of  data-model  misfit  (5)  and  reconstructed 
optical  properties.  The  difference  plots  (Figs.  4(b)  and  (d))  demonstrate  that  Eqs.  9  and 
16  are  equal  within  the  limits  of  the  numerical  precision  to  be  expected  (<  10-8  of  the 
L2-norm  value,  after  the  first  few  iterations).  A  similar  analysis  between  the  original  LM 
update  equation  (Eq.  5)  and  its  alternative  form  (Eq.  6)  was  performed  gave  similar  results 
(not  shown  here).  Reconstruction  results  with  different  target  shapes  and  positions  are 
summarized  in  Table- 1  which  reports  the  mean  and  standard  deviation  of  recovered  /.(„  and 
fj!a  values  in  the  background  and  target  areas.  Note  that  the  recovered  optical  properties 
between  Z  =  15  mm  and  Z  =  -15  mm  were  used  because  the  reconstructed  optical  properties 
were  equal  to  the  actual  background  values  and  the  standard  deviation  was  zero  (within 
round-off  error  limits)  above/below  these  Z- values. 

To  show  the  robustness  of  the  GLS  procedure,  data  with  5%  noise  was  used  in  the 
reconstruction  of  a  cylindrical  target  located  in  the  center  (as  shown  first  rows  of  Fig.  5(a) 
and  (b)).  The  reconstruction  results  using  the  LM  and  GLS  schemes  are  presented  in  middle 
and  last  rows  of  Fig.  5(a)  and  (b)  respectively.  The  GLS  minimization  technique  was  able 
to  localize  the  target  more  clearly  than  the  LM  method. 

To  see  the  effect  of  target  size  on  the  recovery  of  its  contrast  using  these  reconstruction 
techniques,  a  series  of  simulations  were  performed  where  the  diameter  of  the  spherical  target 
located  at  the  center  was  varied  (from  10  to  35  mm).  One  set  of  results  is  presented  for  the 
15  mm  diameter  target  in  Fig.  3.  Another  sample  set  for  a  target  with  diameter  of  10  mm 
is  shown  in  Fig.  6.  A  comparison  plot  is  given  in  Fig.  7.  The  data  noise  level  for  the  cases 
considered  here  is  1%.  Increase  in  the  diameter  of  the  spherical  target  increases  the  contrast 
recovery. 
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A  performance  comparison  of  these  reconstruction  techniques  with  increases  in  target 
(spherical  object  with  diameter  of  25  mm)  contrast  (from  2  to  10  with  respect  to  background 
optical  properties)  located  at  the  center  (0,0,0)  and  (20,0,0)  is  presented  in  Fig.  8.  Again, 
the  noise  level  in  the  data  was  1%.  The  recovery  of  contrast  is  much  lower  in  the  case  of 
the  centered  target  compared  to  the  off-centered  location  in  both  LM  and  GLS  techniques. 
Between  the  LM  and  GLS  methods,  the  later  performs  better  in  terms  of  recovery  of  contrast. 

A  study  was  conducted  to  evaluate  estimation  parameter  independence  (cross-talk)  in 
these  reconstruction  procedures.  The  spherical  target  having  a  contrast  of  2:1  only  in  /j,a 
was  considered  at  the  center  and  near  the  boundary  of  the  imaging  domain.  Synthetic  data 
with  1%  noise  was  used  in  the  reconstructions,  and  the  results  are  presented  in  Table-2  in 
terms  of  recovered  mean  and  standard  deviation  of  the  optical  properties  .  The  recovery  of 
contrast  was  higher  in  the  GLS  case  and  the  amount  of  cross-talk  was  less  (roughly  50%  in 
the  LM  case  compared  to  30%  in  the  GLS  case).  A  similar  study  with  a  contrast  of  only 
in  fjfs  was  also  conducted  (not  shown  )and  it  also  showed  a  similar  trend  in  terms  of  cross-talk. 

Finally,  using  the  experimental  data  (Sec.  5.B)  collected  using  a  multi-layered  gelatin 
phantom,  reconstructions  using  both  LM  and  GLS  techniques  were  performed.  For  the  GLS 
technique,  experimental  data  variance  was  estimated  (Sec.  5. A)  using  analytical  equation, 
given  in  Eq.  22.  Two-dimensional  cross-sections  of  the  actual  and  reconstructed  optical  prop¬ 
erties  distribution  are  plotted  in  Fig.  10  (a)  and  (b).  The  middle  and  bottom  rows  corresponds 
to  the  reconstruction  results  obtained  by  LM  and  GLS  techniques.  As  the  variance  in  the 
data  (noise  charterstics)  are  embedded  in  the  GLS  reconstruction  procedure,  resulting  in 
optimal  weighting  (highly  noisy  data  points  gets  less  weightage  and  vice  versa  for  less  noisy 
data  points),  leading  to  better  quantification  of  tumor  region. 

7.  Discussion 

Appendix  (A.l)  presents  a  computationally  efficient  form  for  implementing  an  iterative  GLS 
reconstruction  scheme  which  reduces  the  dimensionality  of  the  matrix  to  be  inverted.  The 
alternative  forms  for  other  minimization  methods  are  also  developed  in  the  appendices  (A. 2 
and  A. 3).  Appendix-A.4  presents  expressions  for  estimating  the  operations  count  of  both 
forms  of  the  GLS  update  equations  (Eqs.  9  and  16)  for  a  single  iteration.  This  appendix 
(A. 4)  also  gives  operations  count  for  both  LM  and  its  alternative  form  (Eqs.  5  and  6)  as 
well.  Fig.  2  shows  a  log-log  plot  of  operations  count  as  a  function  of  the  ratio  of  number 
of  optical  property  parameters  (NN)  to  number  of  measurements  ( NM )  which  determines 
the  form  of  the  GLS  and  LM  update  equation  to  be  preferred  (given  that  both  produce 
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numerically  equivalent  results  as  reported  in  Fig.  4).  For  example,  when  spatial-priors  are 
available,  the  number  of  optical  unknowns  can  be  reduced  to  the  number  of  regions  that 
can  be  segmented32  which  in  the  case  of  breast  tissue,  is  typically  NN  =  3  (assuming  fatty, 
fibroglandular  and  tumor  regions).32  Here,  since  NM  NN,  the  original  LM  and  GLS 
update  equation  (Eqs.  5  and  9  respectively)  is  effective.  In  under-determined  problems,  such 
as  the  cases  considered  in  this  paper,  where  NM  -C  NN  (NN/NM  ratio  of  12  in  the  test 
problems),  the  GLS  alternative  form  reduces  the  number  of  operations  (by  up  to  two  times 
in  Fig.  2(a)  when  NN/NM  =  12).  In  fact,  the  alternative  form  of  the  GLS  update  equation 
(Eq.  16)  becomes  effective  once  NN/NM  >  6  and  the  number  of  operations  decreases  by  an 
order  of  magnitude  when  NN/NM  reaches  100.  For  the  LM  minimization  scheme,  alternative 
form  reduces  the  number  of  operations  by  a  factor  of  6  in  Fig.  2(b)  when  NN/NM  =  12.  The 
memory  required  for  inverting  such  matrices  is  plotted  in  Fig.  2(c)  as  a  function  of  NN/NM. 
It  is  also  important  to  recognize  that  the  memory  required  to  complete  these  operations 
can  become  critical  because  the  cache  sizes  and  RAM  available  on  different  architectures  is 
variable  but  influences  the  efficiency  of  the  computational  processes  executing  on  a  given 
platform. 

Under-determinedness  of  the  imaging  problem  (i.e.  NN/NM  >  1)  leads  to  non-uniqueness 
in  the  solution  space,  but  regularization  helps  to  give  a  unique  solution  in  these  cases. 
Typically  NN/NM  values  are  between  2  to  10  for  a  typical  two-dimensional  (2D)  problem, 
as  the  choice  depends  on  the  expected  resolution  in  the  reconstructed  image,  imaging 
domain  size,  shape,  data-collection  geometry  and  prior  information  available.  For  a  3D 
problem,  this  choice  (NN/NM  values)  also  depends  on  these  factors,  but  one  expects  this 
ratio  to  be  higher  than  the  2D  case  as  the  imaging  domain  size  is  bigger  reflecting  in  the 
number  of  imaging  parameters  (NN)  to  be  larger  (Typical  example:  NN  =  600  (in  2D)  and 
6000  (in  3D)).  Ideally,  one  would  like  to  have  this  ratio  (NN/NM)  constant  between  2D 
and  3D,  which  implies  that  the  number  of  measurements  has  to  increase  by  the  same  factor 
(typical  case  requires  a  factor  of  10),  which  might  not  be  feasible  due  to  instrumentation 
constrains.24  This  leads  to  choice  of  NN/NM  greater  than  6  at  least,  where  the  derived 
forms  are  effective  (even  though  in  the  case  of  LM,  the  alternative  form  is  effective  when 
NN/NM  >  2).  It  is  also  important  to  note  that,  ideally  one  would  like  mesh  the  volume 
where  the  sensitivity  is  greater  (for  the  imaging  domain  discussed  here  z  between  -25  mm 
and  25  mm)  finer  and  rest  of  the  domain  coarser,  to  keep  NN/NM  in  the  same  range  as  in 
2D  (typically  3-8).  This  adaptive  meshing  gets  trickier  when  the  patients  are  imaged,  in 
our  experience,  we  could  not  find  appropriate  tools  that  could  be  used  in  real-time  for  this 
meshing  problem.  Even  though  efforts  to  solve  this  patient-specfic  adaptive  meshing  is  being 
pursued.30,31  Even  when  NN/NM  =  3  (lowest  ratio  one  expects  in  a  3D  imaging  problem), 
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from  Fig.  2(b),  the  alternative  form  (LM-AF,  Eq.  6)  in  the  case  of  LM  minimization  scheme 
becomes  effective. 

Figure  4  demonstrates  that  the  two  GLS  update  forms  (Eqs.  9  and  16)  are  equivalent 
numerically  (within  the  numerical  precision  of  the  L2  norm  value).  For  the  cases  considered 
here,  the  computation  time  for  each  iterative  update  using  Eq.  16  was  approximately  46 
minutes,  which  was  three  times  faster  than  with  Eq.  9  (computation  time  ~  126  minutes 
per  iteration).  In  terms  of  operations  count,  a  factor  of  2  reduction  would  be  estimated 
from  Fig.  2(a)  for  the  NN/NM  ratio  involved.  The  deviation  in  run-time  that  occurs  in 
practice  is  likely  due  to  the  cost  of  memory  management  alluded  to  above  in  the  case  of 
Eq.  9.  It  is  also  important  to  note  that  implementation  of  Eq.  9  requires  an  inversion  of 
the  covariance  matrix  (Eq.  10),  whereas  this  matrix  can  be  used  directly  in  Eq.  16.  In  the 
case  of  Levenberg-Marquardt  (LM)  update  equations,  computation  time  for  an  iteration 
using  Eq.  6  was  approximately  21  minutes  and  for  Eq.  5  was  91  minutes.  The  deviation 
from  factor  of  6  (from  Fig.  2(b))  is  mainly  due  to  the  memory  required  to  perform  these 
operations,  which  affects  run  time  in  turn. 

Middle  and  last  rows  of  Fig.  3(a)  and  (b)  indicate  that  with  1%  noise  in  the  data,  LM  has 
failed  to  recover  /ua  in  a  spherical  target  with  a  diameter  of  15  mm  located  at  the  center, 
whereas  GLS  was  able  to  identify  the  target  very  well.  The  failure  of  LM  minimization  is 
indicative  of  a  lack  of  sensitivity  at  the  center  of  the  domain24  which  is  improved  through 
the  GLS  approach  by  including  the  noise  characteristics  and  covariances  associated  with  the 
problem.  When  the  same  target  is  located  near  the  boundary  (at  (30,0,0)),  both  techniques 
were  able  to  recover  the  contrast  approximately  20%  better  relative  to  the  center  position 
(Table- 1). 

When  1%  noisy  data  generated  from  a  centered  cylindrical  target  with  a  diameter  of  15 
mm  (first  rows  of  Fig.  5(a)  and  (b))  was  used  in  the  reconstruction,  both  LM  and  GLS 
techniques  were  able  to  recover  approximately  50%  of  the  expected  contrast  (Table-1).  For 
the  same  type  of  target  located  near  the  boundary  (at  (30,0)),  the  recovery  of  contrast  was 
approximately  70%  in  the  case  of  GLS.  For  LM,  the  recovery  was  only  50%  for  /ia  and  85% 
for  /%  under  these  conditions.  The  reconstruction  results  also  show  that  recovery  of  the 
centered  target  is  always  poor  relative  to  an  object  near  the  boundary.  This  is  primarily  due 
to  the  hyper-sensitivity  at  the  boundary  in  these  cases.24  The  extended  cylindrical  target 
is  essentially  equivalent  to  the  two-dimensional  case  of  a  circular  inclusion  and  the  trend 
observed  in  3D  of  recovering  more  contrast  for  a  target  near  the  boundary  is  similar  to  the 
behavior  found  in  2D.24 
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When  5%  noisy  data  was  used,  LM  reconstruction  (middle  rows  of  Fig.  5(a)  and  (b)) 
performs  poorly  in  terms  of  localization  of  the  target,  whereas  GLS  was  able  to  reconstruct 
optical  images  with  better  quality  and  quantitation  (up  to  70%).  Even  though  the  recon¬ 
structed  results  using  very  noisy  data  were  presented  here  from  only  one  type  of  target, 
similar  trends  were  also  observed  in  other  cases  that  mimic  the  2D  reconstructions  reported 
in  Ref.17  These  results  show  that  GLS  outperforms  LM  even  though  the  data  noise  level 
is  high  because  stability  is  retained  by  including  the  noise  characteristics  into  the  weight 
matrices  used  for  normalization. 

Accuracy  in  contrast  recovery  of  local  targets  increases  as  the  size  of  the  target  increases, 
as  shown  in  Fig.  7.  For  example,  the  contrast  recovered  for  a  centered  target  below  20 
mm  in  diameter  is  only  about  30%  of  the  true  value  and  as  low  as  10%  for  the  LM 
algorithm  (Fig.  7),  whereas  increasing  the  size  of  target  to  30  mm  leads  to  quantitative 
accuracy  near  100%.  The  GLS  approach  provides  maximal  contrast  recovery  and  superior 
image  quality  at  all  sizes  relative  to  LM  (example:  Fig.  6).  Even  when  the  target  size  is  as 
low  as  5  mm  (Fig.  6),  the  object  was  well  localized  in  the  GLS  case,  but  not  with  LM  (Fig.  6). 

The  performance  comparison  of  the  algorithms  in  terms  of  contrast  recovery  (Fig.  8) 
confirms  that  the  position  of  the  target  dictates  the  response.  When  the  target  had  10:1 
contrast  in  comparison  to  the  background,  the  maximum  recovery  of  contrast  was  ~  5:1. 
GLS  outperformed  LM  in  this  regard  but  there  is  a  plateau  in  recovery  of  contrast  at  400% 
of  the  background  value. 

Table-2  shows  that  estimation  parameter  dependence  (cross-talk)  is  lower  (by  20%)  for 
GLS  compared  to  LM,  by  reinforcing  the  independence  of  ft,a  and  //(  through  elimination  of 
any  cross-terms  in  the  weight  matrix  WM_Mo  (Eq.  15).  The  inter-parameter  dependence  is 
complex  because  of  the  non-unique  relationship  between  the  optical  property  distribution 
and  the  incomplete  boundary  data,  indicating  that  different  formulations  of  the  inversion 
tend  to  perform  differently.  Nonetheless,  the  estimation  parameter  dependence  is  substan¬ 
tially  higher  in  3D  data-limited  situations,  relative  to  when  the  ratio  of  data  to  number  of 
estimates  is  less  skewed. 

In  the  case  of  phantom  data  (Fig.  10),  as  expected  GLS  reconstructions  showed  more 
promising  results  in  this  test  case.  Characterization  of  the  data  collection  system,  leading 
to  variance  estimation  depending  on  the  voltage  measure  at  PMT  (Fig.  9)  enabled  the 
employment  of  the  GLS  technique  for  experimental  data  reconstructions.  Both  techniques 
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(LM  and  GLS)  were  able  to  give  qualitative  information  about  the  target,  in  terms  of 
quantification  the  GLS  technique  overtakes  the  LM  technique  (Fig.  10).  It  should  be  noted 
that  this  type  of  charcterization  of  the  experimental  system  does  not  take  into  account 
coupling  errors  between  the  light  collection/delivery  fiber  and  tissue  surface.  These  kind 
of  unsystematic  errors  are  difficult  to  estimate  as  it  depends  on  many  parameters,  such  as 
tissue  surface  roughness,  tissue  elastic  properties,  design  of  tissue-fiber  coupling  interface, 
repetability,  and  alingment  of  fibers.  But  development  of  these  type  of  methods  including 
systematic  noise  charcterstics  in  the  reconstruction  procedure  will  be  taking  a  step  in  the 
right  direction.  Moreover,  as  it  can  be  seen  from  Eq.  22,  covarinces  among  the  data  points 
was  ignored,  making  the  data  weight  matrix  (W<s)  a  simple  diagonal  matrix.  Inclusion 
of  covariances  can  offer  a  better  weighting  in  case  of  experimental  data,  this  is  under 
investigation  as  of  now.  Even  though  this  test  case  showed  very  promising  results,  in  our 
experience,  in  cases  where  coupling  errors  are  dominant  in  the  data,  the  GLS  scheme  did  not 
yield  any  meaningful  results.  In  these  cases,  the  LM  technique  was  able  to  give  reasonable 
results. 


Even  though  the  data  used  here  is  generated  by  using  in-plane  data,  it  was  collected  only 
from  the  source  fiber  plane,  previous  invistigations  indicate  that  the  use  of  out-of-plane  data 
(when  the  data  was  collected  from  rest  of  the  fibers  in  all  three  planes)  may  not  give  enough 
advantage  in  terms  of  reconstructed  image  quality  given  an  increase  in  the  data-acquistion 
time  and  computational  cost.24  It  is  also  important  to  note  that,  the  results  presented  here 
uses  a  cylindrical  imaging  domain,  this  study  is  generic  in  nature,  especially  in  terms  of 
proving  the  computational  efficiency  of  alternative  forms  (Fig.  2),  as  NN/NM  was  changed 
over  a  range  of  0.0028  to  100  (spanning  from  well-determined  to  highly  under-determined 
problems) . 


Partial  volume  effects  can  be  observed  in  the  recovery  of  contrast  as  a  function  of  target 
size.  The  recovery  of  contrast  was  much  higher  for  the  extended  cylinder  target  compared 
to  the  spherical  inclusion  (Table-1).  The  quantitative  accuracy  of  reconstructed  images  in¬ 
creases  with  an  increase  in  target  size  (Fig.  7).  GLS  reconstruction  results  of  the  data  from  a 
centered  cylindrical  object  are  encouraging,  demonstrating  recovery  of  more  than  30%  con¬ 
trast  in  this  case  (other  Newton-type  algorithms  have  reported  a  maximum  of  10%  contrast 

7  11  1 3  1  4  24  \ 

recovery  5  5  5  5  ). 


8.  Conclusions 

Three-dimensional  diffuse  optical  tomography  is  more  computationally  intensive  because  of 
the  size  of  the  parameter  space  to  be  reconstructed.  Newton-based  inversion  methods  that 
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operate  on  a  Hessian  matrix,  which  has  dimensions  of  the  number  of  measurements  rather 
than  the  number  of  parameters,  can  be  derived  using  the  Sherman-Morrison- Woodbury 
identity  and  become  computationally  more  efficient  once  the  number  of  estimation  parame¬ 
ters  exceeds  two  times  the  number  of  measurements.  Representative  exampies  demonstrate 
that  this  form  of  update  equation  can  be  at  least  six  times  faster  in  practice  in  the 
highly  under-determined  problems  which  commonly  occur  in  3D.  Three  dimensional  diffuse 
optical  tomographic  reconstruction  algorithms  also  suffer  from  partial  volume  effects  that 
degrade  significantly  the  accuracy  with  which  optical  properties  can  be  quantified.  The 
GLS  approach  which  incorporates  structured  weight  matrices  consisting  of  the  variance  and 
covariance  of  the  data-model  misfit  and  the  optical  properties,  improves  the  quantification 
of  optical  properties  by  at  least  20%  in  3D.  The  GLS  estimate  is  also  robust  to  data  noise 
as  high  as  5%  -  conditions  under  which  other  algorithms  fail  when  the  problem  is  highly 
under-determined.  By  characterizing  the  detector  noise  for  systematic  errors,  using  a  multi¬ 
layered  gelatin  phantom  data,  the  GLS  technique  can  be  easily  employed  for  reconstructing 
experimental  data  and  can  yield  better  quantification  of  targets  compared  to  conventional 
reconstruction  methods.  Future  investigations  will  include  thorough  examination  of  the 
GLS  technique  when  applied  to  phantom  and  clinical  data  and  extension  of  the  technique 
to  direct-spectral  reconstruction.  The  test  data  used  in  this  article,  along  with  computer 
algorithms,  are  available  on  a  web  page.33 

A.  Appendix 


A.l.  Alternative  form  for  GLS  update  equation 

Before  deriving  the  alternative  form,  it  is  useful  to  catalog  several  properties  of  the 
weight-matrices  and  their  inverses.28 


Wj^Cj)-1;  W^M  =  (tv*,)-1 
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(n,\T  —  m.  in  \t  _  n 


If  a  square  matrix  A,  has  block  form 
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with  dimensions  (2NM+2NN)x(2NM+2NN),  it  is  readily  shown  to  be  symmetric  by  invoking 
the  relationships  in  Eq.  23: 
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Since  inverses  of  both  C$  and  exist,  then  A  1  also  exists  and  can  be  expressed  in 

block  form  as  well 


in  which  case 
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requires  that 

— QP  +  JR  =  I 
— QQ  +  JS  =  0 
JTP  +  \V„  l{  =  0 
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These  relationships  can  be  manipulated  through  a  series  of  substitutions  to  express  the  blocks 
of  A-1  in  terms  of  combinations  of  the  block  components  of  A.  Specifically,  Eqs.  29  and  30 
along  with  the  weight  matrix  properties  in  Eqs.  23  imply  that 


Q  =  IT^JS 


(32) 


and 

R  =  -C^0JTP. 

Substituting  Eq.  33  into  Eq.  28  to  form  the  expression 

P  =  -(C(5  +  JC/i_/i0JT)“1 

which  is  put  back  into  Eq.  33  produces 

r  =  c„jt  (c5  +  jc„jt)_1 

A  similar  series  of  steps  starting  with  Eqs.  32  and  31  to  write 

s  =  (wm_w  +  jtw5j)_1 


which  is  combined  again  with  Eq.  32  yields 

q  =  w5j  (w^  +  jWvj)-1. 

Since  A  is  symmetric  and  invertible,  A-1  is  symmetric  as  well. 
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in  which  case 


QT  =  R  (39) 

Substituting  the  forms  of  Q  and  R  (Eqs.  37  and  35,  respectively)  into  Eq.  39  results  in 


1  T 


W^J  (Wn  +  JTW5J)  1  =  CV^ JT  (Cs  +  3Cn JT) 


T\- 1 


(40) 


Equation  38  also  requires  ST  =  S,  where  S  is  given  by  Eq.  36,  which  when  identified  in  the 
term  on  the  left  side  of  Eq.  40  allows  it  to  be  rewritten  as 
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RJ  +  SW^  =  I. 


Solving  for  S 


and  substituting  Eqs.  35  and  36  for  R  and  S  in  Eq.  44  produces 

(WM_M0  +  JtW5J)_1  =  CV^  -  Cm_wJt  (Cs  +  JCM_MoJ T)_1  JCM_ 
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Note  that  this  derivation  was  adapted  from  Liebelt  et.  al.34  A  variant  of  Eq.  45  exists  in  the 
literature  with  many  names,  the  most-common  being  the  Sherman-Morrison- Woodbury 
identity. 35-42  It  is  also  known  as  the  matrix  inversion  lemma.43,44  Even  though  one 
can  start  from  this  equation  and  derive  the  alternative  forms,  the  complete  derivation  is 
presented  here  for  completeness. 

Substituting  Eq.  45  back  into  Eq.  9  yields 
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,  -l 


(Mi  — 1  Mo)} 

(46) 

3  (Mi-i  —  Mo)} 

(47) 

-i  -  Mo)} 

(48) 

as  the  alternative  form  for  the  update  equation  (Eq.  16) 


20 


The  next  two  subsections  show  the  alternative  forms  of  other  least-squares  minimization 
techniques,  namely  Levenberg-Marquardt  and  Tikhonov  minimizations. 


A. 2.  Alternative  form  for  LM  update  equation 

The  Levenberg-Marquardt  (LM)  update  equation  (Eq.  5)  is  a  special  case  of  the  GLS 
update  equation  (Eq.  9)  when  WM_M0  =  al  and  W §  =  I  (see  Sec.  III.B.4  in  Ref.17).  Using 
these  forms  in  Eq.  41  leads  to  an  alternative  form  to  Eq.  5 

A im  =  (crI)-1JT  (J(aI)"1JT  +  I"1)-1  i  (49) 


Rearranging  the  terms  in  Eq.  49  leads  to 


-l 


a  \  a  ) 

which  can  be  simplified  to  produce 

A m  =  JT  (JJT  +  o;l)_1  8i-\. 

Eq.  51  is  also  known  as  under-determined  form  in  the  literature.11,14 
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A. 3.  Alternative  form  for  Tikhonov  update  equation 
The  objective  function17,45  for  the  Tikhonov  scheme  is 

l|y-G'(/r)||2  +  A||L(/i-p0)||2  (52) 


Minimization  of  Eq.  52  and  linearizing  the  problem  leads  to  update  equation17,32 

A/q  =  [JTJ  +  ALrL] _1  (J TSi.1  -  ALTL(/q_!  -  p0))  .  (53) 

Equation  53  is  a  special  case  of  the  GLS  update  equation  (Eq.  9)  with  weight  matrices  (see 
Sec.  III.B.4  in  Ref.17) 

Wi  =  I;  W„_w  =  ALtL  (54) 


From  Eq.  48  one  can  write 


A  pi  = 


I  —  (ALrL)  1Jt(^J(ALtL)  1Jt  +  I“1)  J  {(ALtL)  1  JTI<Si_i  -  (W_1  -Po)}  (55) 


which  leads  to 


A  m  = 


I-(LtL)  1Jt(j(LtL)  XJt  +  Ai)  {(ALtL)  1  JTSi-1  -  (/Zi_i  -Mo)}  (56) 


Assuming  /q_i  =  //0 ,  the  single-step  Tikhonov  update  equation  (or  its  equivalent)14,46,47 
becomes 

A  m  =  (JT  J  +  ALtL)  _1  J  T8i-l  (57) 
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-1 


(58) 


Using  Eqs.  41  and  54  leads  to 


A  &  =  (ALtL)  1  JT  (j  (ALtL)  1  JT  +  I"1) 


5 


l-l 


which  can  be  rearranged  to 

A  pi  =  (LtL)  _1  JT  (j  (LtL)  _1  JT  +  Al)  _1  5^!  (59) 

Eq.  59  is  also  known  as  under-determined  Tikhonov  single-step  update  equation.11, 14,16 


A. 4-  Calculation  of  number  of  operations  for  LM  and  GLS  update  equations 

The  number  of  operations  was  estimated  by  assuming  that  divisions/multiplications 
consume  most  of  the  processor  cycles.  Note  that  Gaussian  elimination  was  used  in 
computing  matrix  inversion.  Typically,  Gaussian  elimination  for  an  NxN  matrix  requires 
((iV3/ 3)  +  N 2  —  (N/ 3))  operations.42  The  memory  required  to  invert  a  matrix  of  dimension 
NxN  is  N2.42  The  number  of  operations  only  includes  solution  of  the  update  equation  and 
does  not  account  for  the  number  of  operations  required  to  form  the  matrices/ vectors  used 
in  these  equations. 


For  the  GLS  update  equation  (Eq.  9),  the  number  of  operations  required  for  iteration  i  is 
(using  the  dimensions  defined  after  Eq.  15) 

Number  of  operations  =  [(2 NM  *  2 NM  *  2NN)  +  (2 NN  *  2 NM  *  2 NN)  +  ( +  (2 NN)2  -  ^^)]  ^ 

+(2NN  *  2NN  *  1)  +  [(27V7V  *  2NM  *  1)  +  (2 NM  *  2NM  *  1)  +  (27V TV  *  2NN  *  1)] 

For  the  alternative  form  for  the  GLS  update  equation  (Eq.  16),  the  number  of  operations  is 


Number  of  operations  =  [(27V  AT  *  2NN  *  2 NM)  +  (2 NM  *  2NN  *  2 NM)  +  (27V7V  *  2NN  *  2NM) 

+  ( (2jv3M>3  +  (27V AT)2  -  27|M)  +  (2NM  *  2NM  *  2 TV TV)]  (61) 

+(27V7V  *  2 TV TV  *  1)  +  [(27V TV  *  27V TV  *  27VM)  +  (27V TV  *  27 VM  *  1)  +  (27 VM  *  27 VM  *  1)] 


Similarly,  for  the  LM  update  equation  (Eq.  5),  the  number  of  operations  required  for  iteration 
i  is 


Number  of  operations  = 


( (2NN)3  o  2NN\~\ 

(27V AT  *  2NM  *  27V AT)  +  ( 1 +  (27V AT)2 - —  J 


■  (27V AT  *  2NM  *  1) 


(62) 


The  number  of  operations  for  the  alternative  form  for  the  LM  update  equation  (Eq.  6)  is 

Number  of  operations  =  (27V AT  *  2NM  *  1)  +  [(2AT7W  *  27V7V  *  2NM)  +  ( (2jvgM)3  +  ( 2NM )2  - 

+(2NM  *  2NM  *  1) 


(63) 


Note  that  the  computation  time  for  these  update  equations  is  not  only  dependent  on 
the  number  of  operations  needed  to  be  performed,  but  also  on  the  memory  required  for 
implementing  the  operations. 
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Figure  10:  Actual  and  reconstructed  (a)  /ia  and  (b)  //(  distributions  of  a  cylindrical 
target  using  experimental  multi-layered  phantom  data.  Two-dimensional  cross-sections  of 
the  3D  volume  in  2.5  mm  increments  spanning  from  z  =  -12.5  mm  to  z  =  12.5  mm  (from  left 
to  right)  are  shown.  Actual  distributions  are  given  in  the  first  row.  Reconstructed  distribu¬ 
tion  using  the  Levenberg-Marquardt  (LM)  minimization  scheme  and  GLS  minimization 
scheme  are  presented  in  the  middle  and  last  rows  respectively. 
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Tables 


Methods 

Target 

Background 

Target 

shape 

position 

ha 

/4 

ha 

Ms 

Actual 

- 

- 

0.01 

1.0 

0.02 

2.0 

LM 

Spherical 

(0,0,0) 

0.0101±0.0003 

1.0079±0.0322 

0.0104±0.0002 

1.1259±0.0160 

(30,0,0) 

0.0101±0.0006 

1.0063±0.0500 

0.0126±0.0009 

1.4514±0.1152 

Cylindrical 

(0,0) 

0.0102±0.0010 

1.0120±0.0874 

0.0151±0.0012 

1.4308±0.0854 

(30,0) 

0.01014=0.0006 

1.0030±0.0663 

0.0148±0.0015 

1.8406±0.2470 

GLS 

Spherical 

(0,0,0) 

0.01014=0.0001 

1.0100±0.0212 

0.0122±0.0004 

1.2903±0.0326 

(30,0,0) 

0.0100±0.0004 

1.0108±0.0250 

0.0141±0.0006 

1.4498±0.0853 

Cylindrical 

(0,0) 

0.0102±0.0008 

1.00554=0.0500 

0.0159±0.0009 

1.4750±0.0941 

(30,0) 

0.01014=0.0008 

1.0043±0.0588 

0.0170±0.0012 

1.6793=b0.1848 

Table  1 
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Methods 

Target 

position 

Background 

Target 

Ha 

Ms 

Ha 

Ms 

Actual 

- 

0.01 

1.0 

0.02 

1.0 

LM 

(0,0,0) 

0.0101±0.0003 

1.0025±0.0211 

0.0109±0.0001 

1.05004=0.0071 

(30,0,0) 

0.0101±0.0004 

1.0009=b0.0205 

0.0119±0.0004 

1.09344=0.0269 

GLS 

(0,0,0) 

0.0101±0.0001 

1.0016=b0.0269 

0.0126±0.0002 

1.09244=0.0253 

(30,0,0) 

0.0100±0.0004 

1.0029=b0.0189 

0.0136±0.0003 

1.10024=0.0294 

Table  2 
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Figure  1 
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Figure  2 
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Figure  3(a) 
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Figure  3(b) 
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Figure  4 
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GLS  p  LM  P  Actual 
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Figure  5(a) 
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Figure  5(b) 
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GLS  LM  Actual 
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Figure  6(a) 
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GLS  p  LM  Actual 
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Figure  6(b) 
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u 

'a 


Spherical  target  diameter  (in  mm) 


Spherical  target  diameter  (in  mm) 


Figure  7 
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Figure  8 
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o(0)  in  degrees  a(V)/V  (in  %) 


Measured  PMT  voltage  (V,  in  volts)  in  log  scale 


Measured  PMT  voltage  (V,  in  volts)  in  log  scale 

Figure  9 


42 


GLS  LM  Actual 
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Figure  10(a) 
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Figure  10(b) 
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